From 047e805bcf8fa89aadf76bab8e46747fec0f2a80 Mon Sep 17 00:00:00 2001
From: Joris van Zwieten <zwieten@astron.nl>
Date: Mon, 30 Jul 2012 13:11:37 +0000
Subject: [PATCH] Task #3223: Cleaned up Demixer.cc, removed boost::multi_array
dependency, several other source files cleaned up and documented.
---
.gitattributes | 2 +
CEP/DP3/DPPP/include/DPPP/BBSExpr.h | 111 ---
CEP/DP3/DPPP/include/DPPP/CMakeLists.txt | 2 +-
CEP/DP3/DPPP/include/DPPP/Demixer.h | 59 +-
CEP/DP3/DPPP/include/DPPP/EstimateMixed.h | 2 +-
.../DPPP/include/DPPP/ModelComponentVisitor.h | 2 -
CEP/DP3/DPPP/include/DPPP/Patch.h | 18 +-
CEP/DP3/DPPP/include/DPPP/PointSource.h | 9 +-
CEP/DP3/DPPP/include/DPPP/Simulator.h | 75 ++
CEP/DP3/DPPP/src/Apply.cc | 85 +-
CEP/DP3/DPPP/src/BBSExpr.cc | 221 -----
CEP/DP3/DPPP/src/CMakeLists.txt | 2 +-
CEP/DP3/DPPP/src/Demixer.cc | 804 ++++++++----------
CEP/DP3/DPPP/src/EstimateMixed.cc | 324 ++++---
CEP/DP3/DPPP/src/Patch.cc | 43 +-
CEP/DP3/DPPP/src/PhaseShift.cc | 48 +-
CEP/DP3/DPPP/src/PointSource.cc | 21 +-
CEP/DP3/DPPP/src/Simulate.cc | 475 ++++-------
CEP/DP3/DPPP/src/Simulator.cc | 304 +++++++
CEP/DP3/DPPP/src/SourceDBUtil.cc | 5 +-
CEP/DP3/DPPP/src/SubtractMixed.cc | 5 +-
21 files changed, 1160 insertions(+), 1457 deletions(-)
delete mode 100644 CEP/DP3/DPPP/include/DPPP/BBSExpr.h
create mode 100644 CEP/DP3/DPPP/include/DPPP/Simulator.h
delete mode 100644 CEP/DP3/DPPP/src/BBSExpr.cc
create mode 100644 CEP/DP3/DPPP/src/Simulator.cc
diff --git a/.gitattributes b/.gitattributes
index 38185193884..26ba1a72404 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -736,6 +736,7 @@ CEP/DP3/DPPP/include/DPPP/PointSource.h -text
CEP/DP3/DPPP/include/DPPP/Position.h -text
CEP/DP3/DPPP/include/DPPP/RunDetails.h -text
CEP/DP3/DPPP/include/DPPP/Simulate.h -text
+CEP/DP3/DPPP/include/DPPP/Simulator.h -text
CEP/DP3/DPPP/include/DPPP/SourceDBUtil.h -text
CEP/DP3/DPPP/include/DPPP/Stokes.h -text
CEP/DP3/DPPP/include/DPPP/SubtractMixed.h -text
@@ -763,6 +764,7 @@ CEP/DP3/DPPP/src/PointSource.cc -text
CEP/DP3/DPPP/src/Position.cc -text
CEP/DP3/DPPP/src/RunDetails.cc -text
CEP/DP3/DPPP/src/Simulate.cc -text
+CEP/DP3/DPPP/src/Simulator.cc -text
CEP/DP3/DPPP/src/SourceDBUtil.cc -text
CEP/DP3/DPPP/src/Stokes.cc -text
CEP/DP3/DPPP/src/SubtractMixed.cc -text
diff --git a/CEP/DP3/DPPP/include/DPPP/BBSExpr.h b/CEP/DP3/DPPP/include/DPPP/BBSExpr.h
deleted file mode 100644
index 756c17114ed..00000000000
--- a/CEP/DP3/DPPP/include/DPPP/BBSExpr.h
+++ /dev/null
@@ -1,111 +0,0 @@
-//# BBSExpr.h: Create the expression tree for BBS
-//#
-//# Copyright (C) 2012
-//# ASTRON (Netherlands Institute for Radio Astronomy)
-//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
-//#
-//# This file is part of the LOFAR software suite.
-//# The LOFAR software suite 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 3 of the License, or
-//# (at your option) any later version.
-//#
-//# The LOFAR software suite 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 the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
-//#
-//# $Id$
-
-#ifndef DPPP_BBSEXPR_H
-#define DPPP_BBSEXPR_H
-
-// \file
-// Predict visibilities given a source model
-
-#include <DPPP/DPBuffer.h>
-#include <DPPP/DPInfo.h>
-#include <DPPP/DPInput.h>
-
-#include <BBSKernel/MeasurementExprLOFAR.h>
-#include <BBSKernel/CorrelationMask.h>
-#include <BBSKernel/ParmManager.h>
-#include <BBSKernel/VisBuffer.h>
-#include <BBSKernel/Estimate.h>
-#include <ParmDB/SourceDB.h>
-#include <ParmDB/Grid.h>
-
-#include <Common/lofar_vector.h>
-#include <Common/lofar_string.h>
-
-namespace LOFAR {
- namespace DPPP {
-
- // \addtogroup NDPPP
- // @{
-
- class BBSExpr
- {
- public:
- // Construct the expression for the given sky and instrument
- // parameter databases.
- BBSExpr (const DPInfo& info, const string& skyName,
- const string& instrumentName, double elevationCutoff);
-
- ~BBSExpr();
-
- // Set the options.
- void setOptions (const BBS::SolverOptions&);
-
- // Add a model expression for the given source and phase reference.
- void addModel (const string& source, const casa::MDirection& phaseRef);
-
- // Estimate the model parameters.
- void estimate (vector<vector<DPBuffer> >& buffers,
- const BBS::Grid& visGrid, const BBS::Grid& solveGrid,
- const vector<casa::Array<casa::DComplex> >& factors);
-
- // Subtract the sources.
- void subtract (vector<DPBuffer>& buffer, const BBS::Grid& visGrid,
- const vector<casa::Array<casa::DComplex> >& factors,
- uint target, uint nSources);
-
- private:
- // For now, forbid copy construction and assignment.
- BBSExpr(const BBSExpr& other);
- BBSExpr& operator= (const BBSExpr& other);
-
- // Clear the solvables in the model expressions.
- void clearSolvables();
-
- // Set the solvables in the model expressions to the gains.
- void setSolvables();
-
- // Gather information about the instrument from the input meta-data.
- void initInstrument(const DPInfo& info);
-
- //# Data members
- boost::shared_ptr<BBS::SourceDB> itsSourceDB;
- vector<BBS::MeasurementExpr::Ptr> itsModels;
- BBS::Instrument::Ptr itsInstrument;
- casa::MDirection itsDelayRef;
- casa::MDirection itsTileRef;
- BBS::BaselineSeq itsBaselines;
- BBS::CorrelationSeq itsCorrelations;
- BBS::BaselineMask itsBaselineMask;
- BBS::CorrelationMask itsCorrelationMask;
- double itsElevationCutoff;
- BBS::EstimateOptions itsOptions;
- BBS::ParmGroup itsParms;
- vector<BBS::ParmGroup> itsModelParms;
- };
-
-// @}
-
- } //# namespace DPPP
-} //# namespace LOFAR
-
-#endif
diff --git a/CEP/DP3/DPPP/include/DPPP/CMakeLists.txt b/CEP/DP3/DPPP/include/DPPP/CMakeLists.txt
index c46e3612d82..5ecc418523b 100644
--- a/CEP/DP3/DPPP/include/DPPP/CMakeLists.txt
+++ b/CEP/DP3/DPPP/include/DPPP/CMakeLists.txt
@@ -10,7 +10,7 @@ set(inst_HEADERS
StationAdder.h Filter.h
PhaseShift.h Demixer.h
Cursor.h CursorUtilCasa.h Position.h Stokes.h SourceDBUtil.h
- Apply.h EstimateMixed.h Simulate.h SubtractMixed.h Baseline.h
+ Apply.h EstimateMixed.h Simulate.h Simulator.h SubtractMixed.h Baseline.h
ModelComponent.h PointSource.h GaussianSource.h Patch.h
ModelComponentVisitor.h
FlaggerStatistics.h)
diff --git a/CEP/DP3/DPPP/include/DPPP/Demixer.h b/CEP/DP3/DPPP/include/DPPP/Demixer.h
index 5be55300e1d..d001f699b25 100644
--- a/CEP/DP3/DPPP/include/DPPP/Demixer.h
+++ b/CEP/DP3/DPPP/include/DPPP/Demixer.h
@@ -69,8 +69,6 @@ namespace LOFAR {
// Parameters are obtained from the parset using the given prefix.
Demixer (DPInput*, const ParSet&, const string& prefix);
- virtual ~Demixer();
-
// Process the data.
// It keeps the data.
// When processed, it invokes the process function of the next step.
@@ -92,11 +90,6 @@ namespace LOFAR {
virtual void showTimings (std::ostream&, double duration) const;
private:
- void initUnknowns();
-
- // Solve gains and subtract sources.
- void demix();
-
// Add the decorrelation factor contribution for each time slot.
void addFactors (const DPBuffer& newBuf,
casa::Array<casa::DComplex>& factorBuf);
@@ -114,18 +107,8 @@ namespace LOFAR {
vector<MultiResultStep*> avgResults,
uint resultIndex);
- // Calculate the P matrix telling how to deal with sources that will
- // not be predicted.
- // Those sources are the last columns in the demixing matrix.
- vector<casa::Array<casa::DComplex> > getP
- (const vector<casa::Array<casa::DComplex> >& factors, uint nsources);
-
- // Convert a double value to a string (with sufficient precision).
- string toString (double value) const;
-
- // Convert a angle string with an optional unit to radians.
- // The default input unit is degrees.
- double getAngle (const casa::String& value) const;
+ // Solve gains and subtract sources.
+ void demix();
// Export the solutions to a ParmDB.
void dumpSolutions();
@@ -149,22 +132,24 @@ namespace LOFAR {
vector<string> itsModelSources;
vector<string> itsExtraSources;
vector<string> itsAllSources;
-// vector<uint> itsCutOffs;
+// vector<double> itsCutOffs;
uint itsNDir;
uint itsNModel;
+ uint itsNStation;
uint itsNBl;
uint itsNCorr;
uint itsNChanIn;
uint itsNTimeIn;
- uint itsNChanOutSubtr;
+ uint itsNTimeDemix;
uint itsNChanAvgSubtr;
uint itsNTimeAvgSubtr;
- uint itsNTimeChunkSubtr;
+ uint itsNChanOutSubtr;
uint itsNTimeOutSubtr;
- uint itsNChanOut;
+ uint itsNTimeChunk;
+ uint itsNTimeChunkSubtr;
uint itsNChanAvg;
uint itsNTimeAvg;
- uint itsNTimeChunk;
+ uint itsNChanOut;
uint itsNTimeOut;
double itsTimeIntervalAvg;
@@ -186,26 +171,22 @@ namespace LOFAR {
//# shape #directions x #directions.
vector<casa::Array<casa::DComplex> > itsFactorsSubtr;
+ PatchList itsPatchList;
+ Position itsPhaseRef;
+ vector<Baseline> itsBaselines;
+ casa::Vector<double> itsFreqDemix;
+ casa::Vector<double> itsFreqSubtr;
+ vector<double> itsUnknowns;
+ vector<double> itsLastKnowns;
+ uint itsTimeIndex;
+ uint itsNConverged;
+
//# Timers.
NSTimer itsTimer;
NSTimer itsTimerPhaseShift;
NSTimer itsTimerDemix;
NSTimer itsTimerSolve;
-
- uint itsNConverged;
- uint itsTimeCount;
- PatchList itsPatchList;
- vector<Baseline> itsBaselines;
- casa::Vector<double> itsFreqDemix;
- casa::Vector<double> itsFreqSubtr;
- uint itsNTimeDemix;
- uint itsNStation;
- Position itsPhaseRef;
- casa::Array<double> itsUnknowns;
-// casa::Array<double> itsErrors;
- casa::Array<double> itsLastKnowns;
-// vector<casa::MeasFrame> itsFrames;
-// vector<casa::MDirection::Convert> itsConverters;
+ NSTimer itsTimerDump;
};
} //# end namespace
diff --git a/CEP/DP3/DPPP/include/DPPP/EstimateMixed.h b/CEP/DP3/DPPP/include/DPPP/EstimateMixed.h
index 5a0a85abed0..10105f04341 100644
--- a/CEP/DP3/DPPP/include/DPPP/EstimateMixed.h
+++ b/CEP/DP3/DPPP/include/DPPP/EstimateMixed.h
@@ -84,7 +84,7 @@ bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
size_t nChannel, const_cursor<Baseline> baselines,
vector<const_cursor<fcomplex> > data, vector<const_cursor<dcomplex> > model,
const_cursor<bool> flag, const_cursor<float> weight,
- const_cursor<dcomplex> mix, double *unknowns, double *errors);
+ const_cursor<dcomplex> mix, double *unknowns);
// @}
} //# namespace DPPP
diff --git a/CEP/DP3/DPPP/include/DPPP/ModelComponentVisitor.h b/CEP/DP3/DPPP/include/DPPP/ModelComponentVisitor.h
index c09aab84973..c45857cd0f9 100644
--- a/CEP/DP3/DPPP/include/DPPP/ModelComponentVisitor.h
+++ b/CEP/DP3/DPPP/include/DPPP/ModelComponentVisitor.h
@@ -34,7 +34,6 @@ namespace DPPP
class PointSource;
class GaussianSource;
-class Patch;
// \addtogroup NDPPP
// @{
@@ -46,7 +45,6 @@ public:
virtual void visit(const PointSource&) = 0;
virtual void visit(const GaussianSource&) = 0;
- virtual void visit(const Patch&) = 0;
};
// @}
diff --git a/CEP/DP3/DPPP/include/DPPP/Patch.h b/CEP/DP3/DPPP/include/DPPP/Patch.h
index b77a0611444..a1494be23c7 100644
--- a/CEP/DP3/DPPP/include/DPPP/Patch.h
+++ b/CEP/DP3/DPPP/include/DPPP/Patch.h
@@ -41,7 +41,7 @@ namespace DPPP
// \addtogroup NDPPP
// @{
-class Patch: public ModelComponent
+class Patch
{
public:
typedef shared_ptr<Patch> Ptr;
@@ -52,11 +52,12 @@ public:
const string &name() const;
virtual const Position &position() const;
- virtual void accept(ModelComponentVisitor &visitor) const;
+
+ size_t nComponents() const;
+ ModelComponent::ConstPtr component(size_t i) const;
private:
typedef vector<ModelComponent::ConstPtr>::const_iterator const_iterator;
-
const_iterator begin() const;
const_iterator end() const;
@@ -92,6 +93,17 @@ inline const Position &Patch::position() const
return itsPosition;
}
+inline size_t Patch::nComponents() const
+{
+ return itsComponents.size();
+}
+
+inline ModelComponent::ConstPtr Patch::component(size_t i) const
+{
+ return itsComponents[i];
+}
+
+
} //# namespace DPPP
} //# namespace LOFAR
diff --git a/CEP/DP3/DPPP/include/DPPP/PointSource.h b/CEP/DP3/DPPP/include/DPPP/PointSource.h
index 00e1ea2f755..f3e9fe0e60e 100644
--- a/CEP/DP3/DPPP/include/DPPP/PointSource.h
+++ b/CEP/DP3/DPPP/include/DPPP/PointSource.h
@@ -38,8 +38,6 @@ namespace LOFAR
namespace DPPP
{
-class ModelComponentVisitor;
-
// \addtogroup NDPPP
// @{
@@ -60,11 +58,7 @@ public:
template <typename T>
void setSpectralIndex(double refFreq, T first, T last);
- void setPolarizedFraction(double fraction);
-
- void setPolarizationAngle(double angle);
-
- void setRotationMeasure(double rm);
+ void setRotationMeasure(double fraction, double angle, double rm);
Stokes stokes(double freq) const;
@@ -81,6 +75,7 @@ private:
double itsPolarizedFraction;
double itsPolarizationAngle;
double itsRotationMeasure;
+ bool itsHasRotationMeasure;
};
// @}
diff --git a/CEP/DP3/DPPP/include/DPPP/Simulator.h b/CEP/DP3/DPPP/include/DPPP/Simulator.h
new file mode 100644
index 00000000000..b4a5b2aeaea
--- /dev/null
+++ b/CEP/DP3/DPPP/include/DPPP/Simulator.h
@@ -0,0 +1,75 @@
+//# Simulator.h: Compute visibilities for different model components types
+//# (implementation of ModelComponentVisitor).
+//#
+//# Copyright (C) 2012
+//# ASTRON (Netherlands Institute for Radio Astronomy)
+//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//#
+//# This file is part of the LOFAR software suite.
+//# The LOFAR software suite 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 3 of the License, or
+//# (at your option) any later version.
+//#
+//# The LOFAR software suite 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 the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//#
+//# $Id$
+
+#ifndef DPPP_SIMULATOR_H
+#define DPPP_SIMULATOR_H
+
+// \file
+// Compute visibilities for different model components types (implementation of
+// ModelComponentVisitor).
+
+#include <DPPP/Baseline.h>
+#include <DPPP/Cursor.h>
+#include <DPPP/ModelComponent.h>
+#include <DPPP/ModelComponentVisitor.h>
+#include <DPPP/Position.h>
+#include <Common/lofar_complex.h>
+#include <Common/lofar_vector.h>
+
+namespace LOFAR
+{
+namespace DPPP
+{
+
+// \addtogroup NDPPP
+// @{
+
+class Simulator: public ModelComponentVisitor
+{
+public:
+ Simulator(const Position &reference, size_t nStation, size_t nBaseline,
+ size_t nChannel, const_cursor<Baseline> baselines,
+ const_cursor<double> freq, const_cursor<double> uvw,
+ cursor<dcomplex> buffer);
+
+ void simulate(const ModelComponent::ConstPtr &component);
+
+private:
+ virtual void visit(const PointSource &component);
+ virtual void visit(const GaussianSource &component);
+
+private:
+ Position itsReference;
+ size_t itsNStation, itsNBaseline, itsNChannel;
+ const_cursor<Baseline> itsBaselines;
+ const_cursor<double> itsFreq, itsUVW;
+ cursor<dcomplex> itsBuffer;
+ vector<dcomplex> itsShiftBuffer, itsSpectrumBuffer;
+};
+
+// @}
+
+} //# namespace DPPP
+} //# namespace LOFAR
+
+#endif
diff --git a/CEP/DP3/DPPP/src/Apply.cc b/CEP/DP3/DPPP/src/Apply.cc
index 761845bd196..cd5dadbfa8a 100644
--- a/CEP/DP3/DPPP/src/Apply.cc
+++ b/CEP/DP3/DPPP/src/Apply.cc
@@ -31,11 +31,6 @@ namespace DPPP
void apply(size_t nBaseline, size_t nChannel, const_cursor<Baseline> baselines,
const_cursor<double> coeff, cursor<dcomplex> data)
{
- dcomplex Jp_00, Jp_01, Jp_10, Jp_11;
- dcomplex Jq_00, Jq_01, Jq_10, Jq_11;
- dcomplex Jq_00_s0, Jq_10_s0, Jq_01_s1, Jq_11_s1, Jq_00_s2, Jq_10_s2,
- Jq_01_s3, Jq_11_s3;
-
for(size_t bl = 0; bl < nBaseline; ++bl)
{
const size_t p = baselines->first;
@@ -43,62 +38,56 @@ void apply(size_t nBaseline, size_t nChannel, const_cursor<Baseline> baselines,
if(p != q)
{
+ // Jones matrix for station P.
coeff.forward(1, p);
- Jp_00 = dcomplex(coeff[0], coeff[1]);
- Jp_01 = dcomplex(coeff[2], coeff[3]);
- Jp_10 = dcomplex(coeff[4], coeff[5]);
- Jp_11 = dcomplex(coeff[6], coeff[7]);
+ const dcomplex Jp_00(coeff[0], coeff[1]);
+ const dcomplex Jp_01(coeff[2], coeff[3]);
+ const dcomplex Jp_10(coeff[4], coeff[5]);
+ const dcomplex Jp_11(coeff[6], coeff[7]);
coeff.backward(1, p);
+ // Jones matrix for station Q, conjugated.
coeff.forward(1, q);
- Jq_00 = dcomplex(coeff[0], -coeff[1]);
- Jq_01 = dcomplex(coeff[2], -coeff[3]);
- Jq_10 = dcomplex(coeff[4], -coeff[5]);
- Jq_11 = dcomplex(coeff[6], -coeff[7]);
+ const dcomplex Jq_00(coeff[0], -coeff[1]);
+ const dcomplex Jq_01(coeff[2], -coeff[3]);
+ const dcomplex Jq_10(coeff[4], -coeff[5]);
+ const dcomplex Jq_11(coeff[6], -coeff[7]);
coeff.backward(1, q);
+ // Compute (Jp x conj(Jq)) * vec(data), where 'x' denotes the
+ // Kronecker product.
for(size_t ch = 0; ch < nChannel; ++ch)
{
- Jq_00_s0 = Jq_00 * (*data);
- Jq_10_s0 = Jq_10 * (*data);
- ++data;
-
- Jq_01_s1 = Jq_01 * (*data);
- Jq_11_s1 = Jq_11 * (*data);
- ++data;
-
- Jq_00_s2 = Jq_00 * (*data);
- Jq_10_s2 = Jq_10 * (*data);
- ++data;
-
- Jq_01_s3 = Jq_01 * (*data);
- Jq_11_s3 = Jq_11 * (*data);
- ++data;
- data -= 4;
-
- *data = Jp_00 * (Jq_00_s0 + Jq_01_s1)
- + Jp_01 * (Jq_00_s2 + Jq_01_s3);
- ++data;
-
- *data = Jp_00 * (Jq_10_s0 + Jq_11_s1)
- + Jp_01 * (Jq_10_s2 + Jq_11_s3);
- ++data;
-
- *data = Jp_10 * (Jq_00_s0 + Jq_01_s1)
- + Jp_11 * (Jq_00_s2 + Jq_01_s3);
- ++data;
-
- *data = Jp_10 * (Jq_10_s0 + Jq_11_s1)
- + Jp_11 * (Jq_10_s2 + Jq_11_s3);
- ++data;
- data -= 4;
-
- // Move to next channel.
+ // Fetch visibilities.
+ const dcomplex xx = data[0];
+ const dcomplex xy = data[1];
+ const dcomplex yx = data[2];
+ const dcomplex yy = data[3];
+
+ // Precompute terms involving conj(Jq) and data. Each term
+ // appears twice in the computation of (Jp x conj(Jq))
+ // * vec(data).
+ const dcomplex Jq_00xx_01xy = Jq_00 * xx + Jq_01 * xy;
+ const dcomplex Jq_00yx_01yy = Jq_00 * yx + Jq_01 * yy;
+ const dcomplex Jq_10xx_11xy = Jq_10 * xx + Jq_11 * xy;
+ const dcomplex Jq_10yx_11yy = Jq_10 * yx + Jq_11 * yy;
+
+ // Compute (Jp x conj(Jq)) * vec(data) from the precomputed
+ // terms.
+ data[0] = Jp_00 * Jq_00xx_01xy + Jp_01 * Jq_00yx_01yy;
+ data[1] = Jp_00 * Jq_10xx_11xy + Jp_01 * Jq_10yx_11yy;
+ data[2] = Jp_10 * Jq_00xx_01xy + Jp_11 * Jq_00yx_01yy;
+ data[3] = Jp_10 * Jq_10xx_11xy + Jp_11 * Jq_10yx_11yy;
+
+ // Move to the next channel.
data.forward(1);
} // Channels.
+
+ // Reset cursor to the beginning of the current baseline.
data.backward(1, nChannel);
}
+ // Move to the next baseline.
data.forward(2);
++baselines;
} // Baselines.
diff --git a/CEP/DP3/DPPP/src/BBSExpr.cc b/CEP/DP3/DPPP/src/BBSExpr.cc
deleted file mode 100644
index 02ff982bd9e..00000000000
--- a/CEP/DP3/DPPP/src/BBSExpr.cc
+++ /dev/null
@@ -1,221 +0,0 @@
-//# BBSExpr.cc: Create the expression tree for BBS
-//#
-//# Copyright (C) 2012
-//# ASTRON (Netherlands Institute for Radio Astronomy)
-//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
-//#
-//# This file is part of the LOFAR software suite.
-//# The LOFAR software suite 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 3 of the License, or
-//# (at your option) any later version.
-//#
-//# The LOFAR software suite 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 the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
-//#
-//# $Id$
-
-#include <lofar_config.h>
-#include <DPPP/BBSExpr.h>
-#include <DPPP/EstimateNDPPP.h>
-
-#include <BBSKernel/MeasurementAIPS.h>
-
-#include <Common/LofarLogger.h>
-#include <Common/lofar_iostream.h>
-#include <Common/lofar_iomanip.h>
-#include <Common/StreamUtil.h>
-
-#include <measures/Measures/MCDirection.h>
-#include <measures/Measures/MCPosition.h>
-#include <measures/Measures/MeasConvert.h>
-
-using namespace casa;
-using namespace LOFAR::BBS;
-
-namespace LOFAR {
- namespace DPPP {
-
- BBSExpr::BBSExpr (const DPInfo& info, const string& skyName,
- const string& instrumentName, double elevationCutoff)
- : itsBaselineMask (true),
- itsCorrelationMask (true),
- itsElevationCutoff (elevationCutoff)
- {
- // Open the sky SourceDB/ParmDB.
- try {
- itsSourceDB = boost::shared_ptr<SourceDB>
- (new SourceDB(ParmDBMeta("casa", skyName)));
- ParmManager::instance().initCategory(SKY, itsSourceDB->getParmDB());
- } catch (Exception& e) {
- THROW(Exception, "Failed to open sky model parameter database: "
- << skyName);
- }
- // Open the instrument ParmDB.
- try {
- ParmManager::instance().initCategory
- (INSTRUMENT, ParmDB(ParmDBMeta("casa", instrumentName)));
- } catch (Exception& e) {
- THROW(Exception, "Failed to open instrument model parameter database: "
- << instrumentName);
- }
- // Create Instrument instance using the information present in DPInfo.
- initInstrument(info);
- // Store reference directions.
- MDirection itsDelayRef = MDirection::Convert(info.delayCenter(),
- MDirection::J2000)();
- MDirection itsTileRef = MDirection::Convert(info.tileBeamDir(),
- MDirection::J2000)();
-
- // Ignore auto-correlations.
- ASSERT(info.getAnt1().size() == info.getAnt2().size());
- for (size_t i=0; i<info.getAnt1().size(); ++i) {
- int ant1 = info.getAnt1()[i];
- int ant2 = info.getAnt2()[i];
- itsBaselines.append(baseline_t(ant1, ant2));
- if (ant1 == ant2) {
- itsBaselineMask.clear (ant1, ant2);
- }
- }
-
- // Use all correlations.
- ASSERT(info.ncorr() == 4);
- itsCorrelations.append(Correlation::XX);
- itsCorrelations.append(Correlation::XY);
- itsCorrelations.append(Correlation::YX);
- itsCorrelations.append(Correlation::YY);
- }
-
- BBSExpr::~BBSExpr()
- {}
-
- void BBSExpr::setOptions (const SolverOptions& lsqOptions)
- {
- // Initialize parameter estimation options.
- itsOptions = EstimateOptions(EstimateOptions::COMPLEX,
- EstimateOptions::L2,
- false,
- 0,
- false,
- ~flag_t(0),
- flag_t(4),
- lsqOptions);
- }
-
- void BBSExpr::addModel (const string& source, const MDirection& phaseRef)
- {
- // NB: The phase reference needs to be taken from the DPInfo object,
- // because it is not necessarily equal to the phase center of the
- // observation: A PhaseShift step may have shifted the phase center
- // somewhere else.
- MDirection phaseRefJ2000 = MDirection::Convert(phaseRef,
- MDirection::J2000)();
- // Define the model configuration.
- ModelConfig config;
- config.setSources (vector<string>(1, source));
- config.setDirectionalGain();
- config.setCache();
- if (itsElevationCutoff > 0) {
- config.setElevationCutConfig (ElevationCutConfig(itsElevationCutoff));
- }
- // TODO: Find a better way to handle the reference frequency. Here we
- // use a bogus reference frequency, which does not matter since we are
- // not using the beam model. But it's still a hack, of course.
- const double refFreq = 60.0e6;
-
- // Create the model expression.
- itsModels.push_back (MeasurementExpr::Ptr
- (new MeasurementExprLOFAR (*itsSourceDB, BufferMap(), config,
- itsInstrument, itsBaselines, refFreq,
- phaseRefJ2000, itsDelayRef, itsTileRef)));
-
- // Determine parameters to estimate and store for later reference.
- vector<string> incl(1, "DirectionalGain:*"), excl;
- ParmGroup parms = ParmManager::instance().makeSubset
- (incl, excl, itsModels.back()->parms());
- itsModelParms.push_back(parms);
-
- // Update list of unique parameters to estimate (models can in principle
- // share parameters).
- itsParms.insert(parms.begin(), parms.end());
- }
-
- void BBSExpr::estimate (vector<vector<DPBuffer> >& buffers,
- const Grid& visGrid, const Grid& solveGrid,
- const vector<Array<DComplex> >& factors)
- {
- // Set parameter domain.
- ParmManager::instance().setDomain(solveGrid.getBoundingBox());
-
- // Force computation of partial derivatives of the model with respect to
- // the parameters of interest.
- setSolvables();
-
- // Estimate parameter values.
- DPPP::estimate (buffers, itsModels, factors, itsBaselines,
- itsCorrelations, itsBaselineMask, itsCorrelationMask,
- visGrid, solveGrid, itsOptions);
-
- // Flush solutions to disk.
- ParmManager::instance().flush();
- }
-
- void BBSExpr::subtract (vector<DPBuffer>& buffer,
- const Grid& visGrid,
- const vector<Array<DComplex> >& factors,
- uint target,
- uint nSources)
- {
- // Ensure computation of partial derivatives is switched off.
- clearSolvables();
-
- // Subtract selected sources from the data.
- LOG_DEBUG_STR("nSources: " << nSources << ", target: " << target);
- DPPP::subtract (buffer, itsModels, factors, itsBaselines,
- itsCorrelations, itsBaselineMask, itsCorrelationMask,
- visGrid, target, nSources);
- }
-
- void BBSExpr::clearSolvables()
- {
- for (uint i=0; i<itsModels.size(); ++i) {
- itsModels[i]->clearSolvables();
- }
- }
-
- void BBSExpr::setSolvables()
- {
- for (uint i=0; i<itsModels.size(); ++i) {
- itsModels[i]->setSolvables (itsModelParms[i]);
- }
- }
-
- void BBSExpr::initInstrument(const DPInfo& info)
- {
- const size_t nStations = info.antennaNames().size();
- vector<Station::Ptr> stations;
- stations.reserve(nStations);
- for (size_t i=0; i<nStations; ++i) {
- // Get station name and ITRF position.
- casa::MPosition position = MPosition::Convert(info.antennaPos()[i],
- MPosition::ITRF)();
-
- // Store station information.
- stations.push_back(Station::Ptr(new Station(info.antennaNames()(i),
- position)));
- }
-
- MPosition position = MPosition::Convert(info.arrayPos(),
- MPosition::ITRF)();
- itsInstrument = Instrument::Ptr(new Instrument("LOFAR", position,
- stations.begin(),
- stations.end()));
- }
-
- } //# namespace DPPP
-} //# namespace LOFAR
diff --git a/CEP/DP3/DPPP/src/CMakeLists.txt b/CEP/DP3/DPPP/src/CMakeLists.txt
index 2ab7135518e..9c3b4f87d6d 100644
--- a/CEP/DP3/DPPP/src/CMakeLists.txt
+++ b/CEP/DP3/DPPP/src/CMakeLists.txt
@@ -13,7 +13,7 @@ lofar_add_library(dppp
StationAdder.cc Filter.cc
PhaseShift.cc Demixer.cc
Position.cc Stokes.cc SourceDBUtil.cc
- Apply.cc EstimateMixed.cc Simulate.cc SubtractMixed.cc
+ Apply.cc EstimateMixed.cc Simulate.cc Simulator.cc SubtractMixed.cc
ModelComponent.cc PointSource.cc GaussianSource.cc Patch.cc
ModelComponentVisitor.cc
FlaggerStatistics.cc
diff --git a/CEP/DP3/DPPP/src/Demixer.cc b/CEP/DP3/DPPP/src/Demixer.cc
index 74e542eede7..236f43c7f1f 100644
--- a/CEP/DP3/DPPP/src/Demixer.cc
+++ b/CEP/DP3/DPPP/src/Demixer.cc
@@ -35,8 +35,6 @@
#include <DPPP/SourceDBUtil.h>
#include <DPPP/SubtractMixed.h>
-#include <DPPP/PointSource.h>
-
#include <ParmDB/Axis.h>
#include <ParmDB/SourceDB.h>
#include <ParmDB/ParmDB.h>
@@ -58,8 +56,6 @@
#include <casa/Arrays/MatrixMath.h>
#include <scimath/Mathematics/MatrixMathLA.h>
-#include <boost/multi_array.hpp>
-
using namespace casa;
namespace LOFAR {
@@ -67,44 +63,52 @@ namespace LOFAR {
using LOFAR::operator<<;
-// namespace {
-// void pack(const vector<size_t> &directions, size_t nBl, size_t nCh,
-// size_t nCr, const_cursor<dcomplex> in, cursor<dcomplex> out);
-// void pack(const vector<size_t> &directions, size_t nStation,
-// const_cursor<double> in, double *out);
-// void unpack(const vector<size_t> &directions, size_t nStation,
-// const double *in, cursor<double> out);
-// Patch::Ptr makePatch(const string &name);
-// } //# end unnamed namespace
+ namespace
+ {
+ string toString (double value);
+ double getAngle (const String& value);
+ } //# end unnamed namespace
Demixer::Demixer (DPInput* input,
const ParSet& parset, const string& prefix)
- : itsInput (input),
- itsName (prefix),
- itsSkyName (parset.getString(prefix+"skymodel", "sky")),
- itsInstrumentName(parset.getString(prefix+"instrumentmodel",
+ : itsInput (input),
+ itsName (prefix),
+ itsSkyName (parset.getString(prefix+"skymodel", "sky")),
+ itsInstrumentName (parset.getString(prefix+"instrumentmodel",
"instrument")),
- itsTargetSource (parset.getString(prefix+"targetsource", string())),
- itsSubtrSources (parset.getStringVector (prefix+"subtractsources")),
- itsModelSources (parset.getStringVector (prefix+"modelsources",
- vector<string>())),
- itsExtraSources (parset.getStringVector (prefix+"othersources",
- vector<string>())),
-// itsCutOffs (parset.getUintVector (prefix+"elevationcutoffs",
-// vector<uint>())),
-/// itsJointSolve (parset.getBool (prefix+"jointsolve", true)),
- itsNTimeIn (0),
- itsNChanAvgSubtr (parset.getUint (prefix+"freqstep", 1)),
- itsNTimeAvgSubtr (parset.getUint (prefix+"timestep", 1)),
- itsNTimeOutSubtr (0),
- itsNChanAvg (parset.getUint (prefix+"demixfreqstep",
- itsNChanAvgSubtr)),
- itsNTimeAvg (parset.getUint (prefix+"demixtimestep",
- itsNTimeAvgSubtr)),
- itsNTimeChunk (parset.getUint (prefix+"ntimechunk", 0)),
- itsNTimeOut (0),
- itsNConverged (0),
- itsTimeCount (0)
+ itsAvgResultSubtr (0),
+ itsTargetSource (parset.getString(prefix+"targetsource", string())),
+ itsSubtrSources (parset.getStringVector (prefix+"subtractsources")),
+ itsModelSources (parset.getStringVector (prefix+"modelsources",
+ vector<string>())),
+ itsExtraSources (parset.getStringVector (prefix+"othersources",
+ vector<string>())),
+// itsCutOffs (parset.getDoubleVector (prefix+"elevationcutoffs",
+// vector<double>())),
+// itsJointSolve (parset.getBool (prefix+"jointsolve", true)),
+ itsNDir (0),
+ itsNModel (0),
+ itsNStation (0),
+ itsNBl (0),
+ itsNCorr (0),
+ itsNChanIn (0),
+ itsNTimeIn (0),
+ itsNTimeDemix (0),
+ itsNChanAvgSubtr (parset.getUint (prefix+"freqstep", 1)),
+ itsNTimeAvgSubtr (parset.getUint (prefix+"timestep", 1)),
+ itsNChanOutSubtr (0),
+ itsNTimeOutSubtr (0),
+ itsNTimeChunk (parset.getUint (prefix+"ntimechunk", 0)),
+ itsNTimeChunkSubtr(0),
+ itsNChanAvg (parset.getUint (prefix+"demixfreqstep",
+ itsNChanAvgSubtr)),
+ itsNTimeAvg (parset.getUint (prefix+"demixtimestep",
+ itsNTimeAvgSubtr)),
+ itsNChanOut (0),
+ itsNTimeOut (0),
+ itsTimeIntervalAvg(0),
+ itsTimeIndex (0),
+ itsNConverged (0)
{
// Get and set solver options.
// itsSolveOpt.maxIter =
@@ -122,18 +126,13 @@ namespace LOFAR {
// itsSolveOpt.useSVD =
// parset.getBool (prefix+"Solve.Options.UseSVD", true);
- // Maybe optionally a parset parameter directions to give the
- // directions of unknown sources.
- // Or make sources a vector of vectors like [name, ra, dec] where
- // ra and dec are optional.
-
ASSERTSTR (!(itsSkyName.empty() || itsInstrumentName.empty()),
"An empty name is given for the sky and/or instrument model");
// Default nr of time chunks is maximum number of threads.
if (itsNTimeChunk == 0) {
itsNTimeChunk = OpenMP::maxThreads();
}
- // Check that time windows fit nicely.
+ // Check that time windows fit integrally.
ASSERTSTR ((itsNTimeChunk * itsNTimeAvg) % itsNTimeAvgSubtr == 0,
"time window should fit final averaging integrally");
itsNTimeChunkSubtr = (itsNTimeChunk * itsNTimeAvg) / itsNTimeAvgSubtr;
@@ -224,56 +223,15 @@ namespace LOFAR {
itsFirstSteps.push_back (targetAvgSubtr);
// while(itsCutOffs.size() < itsNModel) {
-// itsCutOffs.push_back(0);
+// itsCutOffs.push_back(0.0);
// }
// itsCutOffs.resize(itsNModel);
-
-// itsFrames.reserve(OpenMP::maxThreads());
-// itsConverters.reserve(OpenMP::maxThreads());
-// for(size_t i = 0; i < OpenMP::maxThreads(); ++i) {
-// MeasFrame frame;
-// frame.set(input->arrayPos());
-// frame.set(MEpoch());
-// itsFrames.push_back(frame);
-
-// // Using AZEL because AZELGEO produces weird results.
-// itsConverters.push_back(MDirection::Convert(MDirection::Ref(MDirection::J2000),
-// MDirection::Ref(MDirection::AZEL, itsFrames.back())));
-// }
- }
-
- void Demixer::initUnknowns()
- {
- itsUnknowns.resize(IPosition(4, 8, itsNStation, itsNModel,
- itsNTimeDemix));
- itsUnknowns = 0.0;
-
-// itsErrors.resize(IPosition(4, 8, itsNStation, itsNModel,
-// itsNTimeDemix));
-// itsErrors = 0.0;
-
- itsLastKnowns.resize(IPosition(3, 8, itsNStation, itsNModel));
- for(uint dr = 0; dr < itsNModel; ++dr) {
- for(uint st = 0; st < itsNStation; ++st) {
- itsLastKnowns(IPosition(3, 0, st, dr)) = 1.0;
- itsLastKnowns(IPosition(3, 1, st, dr)) = 0.0;
- itsLastKnowns(IPosition(3, 2, st, dr)) = 0.0;
- itsLastKnowns(IPosition(3, 3, st, dr)) = 0.0;
- itsLastKnowns(IPosition(3, 4, st, dr)) = 0.0;
- itsLastKnowns(IPosition(3, 5, st, dr)) = 0.0;
- itsLastKnowns(IPosition(3, 6, st, dr)) = 1.0;
- itsLastKnowns(IPosition(3, 7, st, dr)) = 0.0;
- }
- }
- }
-
- Demixer::~Demixer()
- {
}
void Demixer::updateInfo (const DPInfo& infoIn)
{
info() = infoIn;
+
// Get size info.
itsNStation = infoIn.antennaNames().size();
itsNChanIn = infoIn.nchan();
@@ -309,13 +267,13 @@ namespace LOFAR {
itsNChanAvgSubtr = info().update (itsNChanAvgSubtr, itsNTimeAvgSubtr);
itsNChanOutSubtr = info().nchan();
ASSERTSTR (itsNChanAvg % itsNChanAvgSubtr == 0,
- "Demix averaging " << itsNChanAvg
- << " must be multiple of output averaging "
- << itsNChanAvgSubtr);
+ "Demix averaging " << itsNChanAvg
+ << " must be multiple of output averaging "
+ << itsNChanAvgSubtr);
ASSERTSTR (itsNTimeAvg % itsNTimeAvgSubtr == 0,
- "Demix averaging " << itsNTimeAvg
- << " must be multiple of output averaging "
- << itsNTimeAvgSubtr);
+ "Demix averaging " << itsNTimeAvg
+ << " must be multiple of output averaging "
+ << itsNTimeAvgSubtr);
// Store channel frequencies for the demix and subtract resolutions.
itsFreqDemix = infoDemix.chanFreqs();
itsFreqSubtr = getInfo().chanFreqs();
@@ -327,24 +285,40 @@ namespace LOFAR {
itsPhaseRef = Position(angles.getBaseValue()[0],
angles.getBaseValue()[1]);
- initUnknowns();
+ // Intialize the unknowns.
+ itsUnknowns.resize(itsNTimeDemix * itsNModel * itsNStation * 8);
+ itsLastKnowns.resize(itsNModel * itsNStation * 8);
+ vector<double>::iterator it = itsLastKnowns.begin();
+ vector<double>::iterator it_end = itsLastKnowns.end();
+ while(it != it_end)
+ {
+ *it++ = 1.0;
+ *it++ = 0.0;
+ *it++ = 0.0;
+ *it++ = 0.0;
+ *it++ = 0.0;
+ *it++ = 0.0;
+ *it++ = 1.0;
+ *it++ = 0.0;
+ }
}
void Demixer::show (std::ostream& os) const
{
os << "Demixer " << itsName << std::endl;
- os << " skymodel: " << itsSkyName << std::endl;
- os << " instrumentmodel:" << itsInstrumentName << std::endl;
- os << " targetsource: " << itsTargetSource << std::endl;
- os << " subtractsources:" << itsSubtrSources << std::endl;
- os << " modelsources: " << itsModelSources << std::endl;
- os << " extrasources: " << itsExtraSources << std::endl;
+ os << " skymodel: " << itsSkyName << std::endl;
+ os << " instrumentmodel: " << itsInstrumentName << std::endl;
+ os << " targetsource: " << itsTargetSource << std::endl;
+ os << " subtractsources: " << itsSubtrSources << std::endl;
+ os << " modelsources: " << itsModelSources << std::endl;
+ os << " extrasources: " << itsExtraSources << std::endl;
+// os << " elevationcutoffs: " << itsCutOffs << std::endl;
// os << " jointsolve: " << itsJointSolve << std::endl;
- os << " freqstep: " << itsNChanAvgSubtr << std::endl;
- os << " timestep: " << itsNTimeAvgSubtr << std::endl;
- os << " demixfreqstep: " << itsNChanAvg << std::endl;
- os << " demixtimestep: " << itsNTimeAvg << std::endl;
- os << " ntimechunk: " << itsNTimeChunk << std::endl;
+ os << " freqstep: " << itsNChanAvgSubtr << std::endl;
+ os << " timestep: " << itsNTimeAvgSubtr << std::endl;
+ os << " demixfreqstep: " << itsNChanAvg << std::endl;
+ os << " demixtimestep: " << itsNTimeAvg << std::endl;
+ os << " ntimechunk: " << itsNTimeChunk << std::endl;
// os << " Solve.Options.MaxIter: " << itsSolveOpt.maxIter << endl;
// os << " Solve.Options.EpsValue: " << itsSolveOpt.epsValue << endl;
// os << " Solve.Options.EpsDerivative: " << itsSolveOpt.epsDerivative << endl;
@@ -379,6 +353,9 @@ namespace LOFAR {
os << " ";
FlagCounter::showPerc1 (os, itsTimerSolve.getElapsed(), self);
os << " of it spent in estimating gains and computing residuals" << endl;
+ os << " ";
+ FlagCounter::showPerc1 (os, itsTimerDump.getElapsed(), self);
+ os << " of it spent in writing gain solutions to disk" << endl;
}
bool Demixer::process (const DPBuffer& buf)
@@ -407,8 +384,7 @@ namespace LOFAR {
}
itsTimerPhaseShift.stop();
- // For each itsNTimeAvg times, calculate the
- // phase rotation per direction.
+ // For each itsNTimeAvg times, calculate the phase rotation per direction.
itsTimerDemix.start();
addFactors (newBuf, itsFactorBuf);
if (itsNTimeIn % itsNTimeAvg == 0) {
@@ -418,29 +394,51 @@ namespace LOFAR {
itsNChanAvg);
// Deproject sources without a model.
deproject (itsFactors[itsNTimeOut], itsAvgResults, itsNTimeOut);
- itsFactorBuf = Complex(); // clear summation buffer
+ itsFactorBuf = Complex(); // Clear summation buffer
itsNTimeOut++;
}
- // Subtract is done with different averaging parameters, so
- // calculate the factors for it.
+ // Subtract is done with different averaging parameters, so calculate the
+ // factors for it.
addFactors (newBuf, itsFactorBufSubtr);
if (itsNTimeIn % itsNTimeAvgSubtr == 0) {
makeFactors (itsFactorBufSubtr, itsFactorsSubtr[itsNTimeOutSubtr],
itsAvgResultSubtr->get()[itsNTimeOutSubtr].getWeights(),
itsNChanOutSubtr,
itsNChanAvgSubtr);
- itsFactorBufSubtr = Complex(); // clear summation buffer
+ itsFactorBufSubtr = Complex(); // Clear summation buffer
itsNTimeOutSubtr++;
}
itsTimerDemix.stop();
- // Do BBS solve, etc. when sufficient time slots have been collected.
+ // Estimate gains and subtract source contributions when sufficient time
+ // slots have been collected.
if (itsNTimeOut == itsNTimeChunk) {
- demix();
+ if(itsNModel > 0) {
+ itsTimerSolve.start();
+ demix();
+ itsTimerSolve.stop();
+ }
+
+ // Clear the input buffers.
+ for (size_t i=0; i<itsAvgResults.size(); ++i) {
+ itsAvgResults[i]->clear();
+ }
+
+ // Let the next step process the data.
+ for (uint i=0; i<itsNTimeOutSubtr; ++i) {
+ itsTimer.stop();
+ getNextStep()->process (itsAvgResultSubtr->get()[i]);
+ itsTimer.start();
+ }
+
+ // Clear the output buffer.
+ itsAvgResultSubtr->clear();
+
+ // Reset counters.
itsNTimeIn = 0;
itsNTimeOut = 0;
itsNTimeOutSubtr = 0;
- itsTimeCount += itsNTimeChunk;
+ itsTimeIndex += itsNTimeChunk;
}
itsTimer.stop();
@@ -449,10 +447,10 @@ namespace LOFAR {
void Demixer::finish()
{
+ itsTimer.start();
+
// Process remaining entries.
if (itsNTimeIn > 0) {
- itsTimer.start();
-
// Finish the initial steps (phaseshift and average).
itsTimerPhaseShift.start();
for (int i=0; i<int(itsFirstSteps.size()); ++i) {
@@ -481,13 +479,36 @@ namespace LOFAR {
// Resize lists of mixing factors to the number of valid entries.
itsFactors.resize(itsNTimeOut);
itsFactorsSubtr.resize(itsNTimeOutSubtr);
+
// Demix the source directions.
- demix();
- itsTimer.stop();
+ if(itsNModel > 0) {
+ itsTimerSolve.start();
+ demix();
+ itsTimerSolve.stop();
+ }
+
+ // Clear the input buffers.
+ for (size_t i=0; i<itsAvgResults.size(); ++i) {
+ itsAvgResults[i]->clear();
+ }
+
+ // Let the next step process the data.
+ for (uint i=0; i<itsNTimeOutSubtr; ++i) {
+ itsTimer.stop();
+ getNextStep()->process (itsAvgResultSubtr->get()[i]);
+ itsTimer.start();
+ }
+
+ // Clear the output buffer.
+ itsAvgResultSubtr->clear();
}
// Write solutions to disk in ParmDB format.
+ itsTimerDump.start();
dumpSolutions();
+ itsTimerDump.stop();
+
+ itsTimer.stop();
// Let the next steps finish.
getNextStep()->finish();
@@ -702,356 +723,222 @@ namespace LOFAR {
factors.reference (newFactors);
}
+ namespace {
+ struct ThreadPrivateStorage
+ {
+ vector<double> unknowns;
+ vector<double> uvw;
+ vector<dcomplex> model;
+ vector<dcomplex> model_subtr;
+ size_t count_converged;
+ };
+
+ void initThreadPrivateStorage(ThreadPrivateStorage &storage,
+ size_t nDirection, size_t nStation, size_t nBaseline, size_t nChannel,
+ size_t nChannelSubtr)
+ {
+ storage.unknowns.resize(nDirection * nStation * 8);
+ storage.uvw.resize(nStation * 3);
+ storage.model.resize(nDirection * nBaseline * nChannel * 4);
+ storage.model_subtr.resize(nBaseline * nChannelSubtr * 4);
+ storage.count_converged = 0;
+ }
+ } //# end unnamed namespace
+
void Demixer::demix()
{
- // Only solve and subtract if multiple directions.
- if (itsNDir > 1) {
- // Collect buffers for each direction.
- vector<vector<DPBuffer> > streams;
- for (size_t i=0; i<itsAvgResults.size(); ++i) {
- // Only the phased shifted and averaged data for directions which have
- // an associated model should be used to estimate the directional
- // response.
- if(i < itsNModel) {
- streams.push_back (itsAvgResults[i]->get());
- }
- // Clear the buffers.
- itsAvgResults[i]->clear();
- }
-
- itsTimerSolve.start();
-
- const size_t nTime = streams[0].size();
- const size_t nThread = OpenMP::maxThreads();
- const size_t nDr = itsNModel;
- const size_t nSt = itsNStation;
- const size_t nBl = itsBaselines.size();
- const size_t nCh = itsFreqDemix.size();
- const size_t nCr = 4;
- const size_t timeFactor = itsNTimeAvg / itsNTimeAvgSubtr;
- vector<DPBuffer> &target = itsAvgResultSubtr->get();
-
-// // Thread-private index per direction of the last known "good" solution.
-// boost::multi_array<int, 2> last(boost::extents[nDr][nThread]);
-// fill(&(last[0][0]), &(last[0][0]) + nDr * nThread, -1);
-
- // Thread-private buffer for the unknowns.
- boost::multi_array<double, 2> unknowns(boost::extents[nThread]
- [nDr * nSt * 8]);
-// boost::multi_array<double, 2> errors(boost::extents[nThread]
-// [nDr * nSt * 8]);
+ const size_t nThread = OpenMP::maxThreads();
+ const size_t nTime = itsAvgResults[0]->get().size();
+ const size_t nTimeSubtr = itsAvgResultSubtr->get().size();
+ const size_t multiplier = itsNTimeAvg / itsNTimeAvgSubtr;
+ const size_t nDr = itsNModel;
+ const size_t nDrSubtr = itsSubtrSources.size();
+ const size_t nSt = itsNStation;
+ const size_t nBl = itsBaselines.size();
+ const size_t nCh = itsFreqDemix.size();
+ const size_t nChSubtr = itsFreqSubtr.size();
+ const size_t nCr = 4;
+ const size_t nSamples = nBl * nCh * nCr;
+
+ vector<ThreadPrivateStorage> threadStorage(nThread);
+ for(vector<ThreadPrivateStorage>::iterator it = threadStorage.begin(),
+ end = threadStorage.end(); it != end; ++it)
+ {
+ initThreadPrivateStorage(*it, nDr, nSt, nBl, nCh, nChSubtr);
// Copy solutions from global solution array to thread private solution
// array (solution propagation between chunks).
- for(size_t i = 0; i < nThread; ++i)
- {
- copy(&(itsLastKnowns(IPosition(3, 0, 0, 0))),
- &(itsLastKnowns(IPosition(3, 0, 0, 0))) + nDr * nSt * 8,
- &(unknowns[i][0]));
- }
-
- // Thread-private buffer for station UVW coordinates.
- boost::multi_array<double, 3> uvw(boost::extents[nThread][nSt][3]);
-
- // Thread-private buffer for the simulated visibilities for all
- // directions at the demix resolution.
- boost::multi_array<dcomplex, 5> buffer(boost::extents[nThread][nDr][nBl]
- [nCh][nCr]);
-
- // Thread-private buffer for the simulated visibilities for a single
- // direction at the residual resolution.
- const size_t nChRes = itsFreqSubtr.size();
- boost::multi_array<dcomplex, 4> residual(boost::extents[nThread][nBl]
- [nChRes][nCr]);
+ copy(itsLastKnowns.begin(), itsLastKnowns.end(), it->unknowns.begin());
+ }
- // Thread-private convergence counter.
- vector<size_t> converged(nThread, 0);
+ const_cursor<double> cr_freq = casa_const_cursor(itsFreqDemix);
+ const_cursor<double> cr_freqSubtr = casa_const_cursor(itsFreqSubtr);
+ const_cursor<Baseline> cr_baseline(&(itsBaselines[0]));
#pragma omp parallel for
- for(size_t i = 0; i < nTime; ++i)
+ for(size_t ts = 0; ts < nTime; ++ts)
+ {
+ const size_t thread = OpenMP::threadNum();
+ ThreadPrivateStorage &storage = threadStorage[thread];
+
+ // Simulate.
+ //
+ // Model visibilities for each direction of interest will be computed
+ // and stored.
+ size_t stride_uvw[2] = {1, 3};
+ cursor<double> cr_uvw_split(&(storage.uvw[0]), 2, stride_uvw);
+
+ size_t stride_model[3] = {1, nCr, nCr * nCh};
+ fill(storage.model.begin(), storage.model.end(), dcomplex());
+ for(size_t dr = 0; dr < nDr; ++dr)
{
- const size_t thread = OpenMP::threadNum();
-
- // Compute elevations and determine which sources are visible.
- // TODO: Investigate segfault in MeasRef...
-// vector<size_t> drUp;
-// Quantum<Double> qEpoch(streams[0][i].getTime(), "s");
-// MEpoch mEpoch(qEpoch, MEpoch::UTC);
-// itsFrames[thread].set(mEpoch);
-
-// for(size_t dr = 0; dr < nDr; ++dr)
-// {
-// MVDirection drJ2000(itsPatchList[dr]->position()[0],
-// itsPatchList[dr]->position()[1]);
-// MVDirection mvAzel(itsConverters[thread](drJ2000).getValue());
-// Vector<Double> azel = mvAzel.getAngle("deg").getValue();
-
-// if(azel(1) >= itsCutOffs[dr])
-// {
-// drUp.push_back(dr);
-// last[dr][thread] = i;
-// }
-// }
-
-// const size_t nDrUp = drUp.size();
-// if(nDrUp == 0)
-// {
-// continue;
-// }
-
- // Simulate.
- size_t strides[3] = {1, nCr, nCh * nCr};
- size_t strides_split[2] = {1, 3};
-
- const_cursor<Baseline> cr_baseline(&(itsBaselines[0]));
- const_cursor<double> cr_freq(&(itsFreqDemix[0]));
- const_cursor<double> cr_freqRes(&(itsFreqSubtr[0]));
-
- cursor<double> cr_split(&(uvw[thread][0][0]), 2, strides_split);
-// for(size_t dr = 0; dr < nDrUp; ++dr)
- for(size_t dr = 0; dr < nDr; ++dr)
- {
- fill(&(buffer[thread][dr][0][0][0]),
- &(buffer[thread][dr][0][0][0]) + nBl * nCh * nCr,
- dcomplex(0.0, 0.0));
-
-// const_cursor<double> cr_uvw =
-// casa_const_cursor(streams[drUp[dr]][i].getUVW());
- const_cursor<double> cr_uvw =
- casa_const_cursor(streams[dr][i].getUVW());
- splitUVW(nSt, nBl, cr_baseline, cr_uvw, cr_split);
-
- cursor<dcomplex> cr_model(&(buffer[thread][dr][0][0][0]), 3,
- strides);
-// simulate(itsPatchList[drUp[dr]]->position(),
-// itsPatchList[drUp[dr]], nSt, nBl, nCh, cr_baseline, cr_freq,
-// cr_split, cr_model);
- simulate(itsPatchList[dr]->position(), itsPatchList[dr], nSt,
- nBl, nCh, cr_baseline, cr_freq, cr_split, cr_model);
- }
-
- // Estimate.
-// vector<const_cursor<fcomplex> > cr_data(nDrUp);
-// vector<const_cursor<dcomplex> > cr_model(nDrUp);
- vector<const_cursor<fcomplex> > cr_data(nDr);
- vector<const_cursor<dcomplex> > cr_model(nDr);
-// for(size_t dr = 0; dr < nDrUp; ++dr)
- for(size_t dr = 0; dr < nDr; ++dr)
- {
-// cr_data[dr] = casa_const_cursor(streams[drUp[dr]][i].getData());
- cr_data[dr] = casa_const_cursor(streams[dr][i].getData());
- cr_model[dr] =
- const_cursor<dcomplex>(&(buffer[thread][dr][0][0][0]), 3,
- strides);
- }
-
- const_cursor<bool> cr_flag =
- casa_const_cursor(streams[0][i].getFlags());
- const_cursor<float> cr_weight =
- casa_const_cursor(streams[0][i].getWeights());
-
- bool status = false;
-// if(nDrUp != nDr)
-// {
-// vector<double> unknowns_packed(nDrUp * nSt * 8);
-// vector<double> errors_packed(nDrUp * nSt * 8);
-// vector<dcomplex> mix_packed(nDrUp * nDrUp * nCr * nCh * nBl);
-
-// // pack unknowns, mixing factors
-// size_t strides_mix[5] = {1, nDrUp, nDrUp * nDrUp, nCr * nDrUp
-// * nDrUp, nCh * nCr * nDrUp * nDrUp};
-// pack(drUp, nBl, nCh, nCr, casa_const_cursor(itsFactors[i]),
-// cursor<dcomplex>(&(mix_packed[0]), 5, strides_mix));
-
-// size_t strides_unknowns[3] = {1, 8, nSt * 8};
-// pack(drUp, nSt, const_cursor<double>(&(unknowns[thread][0]), 3,
-// strides_unknowns), &(unknowns_packed[0]));
-
-// const_cursor<dcomplex> cr_mix(&(mix_packed[0]), 5, strides_mix);
-
-// // estimate
-// status = estimate(nDrUp, nSt, nBl, nCh, cr_baseline, cr_data,
-// cr_model, cr_flag, cr_weight, cr_mix, &(unknowns_packed[0]),
-// &(errors_packed[0]));
-
-// // unpack unknowns, errors
-// unpack(drUp, nSt, &(unknowns_packed[0]),
-// cursor<double>(&(unknowns[thread][0]), 3, strides_unknowns));
-// unpack(drUp, nSt, &(errors_packed[0]),
-// cursor<double>(&(errors[thread][0]), 3, strides_unknowns));
-// }
-// else
- {
- const_cursor<dcomplex> cr_mix =
- casa_const_cursor(itsFactors[i]);
-// status = estimate(nDr, nSt, nBl, nCh, cr_baseline, cr_data,
-// cr_model, cr_flag, cr_weight, cr_mix,
-// &(unknowns[thread][0]), &(errors[thread][0]));
- status = estimate(nDr, nSt, nBl, nCh, cr_baseline, cr_data,
- cr_model, cr_flag, cr_weight, cr_mix,
- &(unknowns[thread][0]), 0);
- }
-
- if(status)
- {
- ++converged[thread];
- }
-
- const size_t nTimeResidual = min(timeFactor, target.size() - i
- * timeFactor);
-
- const size_t nDrRes = itsSubtrSources.size();
-
- size_t strides_model[3] = {1, nCr, nChRes * nCr};
- cursor<dcomplex> cr_model_res(&(residual[thread][0][0][0]), 3,
- strides_model);
+ const_cursor<double> cr_uvw =
+ casa_const_cursor(itsAvgResults[dr]->get()[ts].getUVW());
+ splitUVW(nSt, nBl, cr_baseline, cr_uvw, cr_uvw_split);
+
+ cursor<dcomplex> cr_model(&(storage.model[dr * nSamples]), 3,
+ stride_model);
+ simulate(itsPatchList[dr]->position(), itsPatchList[dr], nSt,
+ nBl, nCh, cr_baseline, cr_freq, cr_uvw_split, cr_model);
+ }
- for(size_t j = 0; j < nTimeResidual; ++j)
- {
-// for(size_t dr = 0; dr < nDrUp; ++dr)
- for(size_t dr = 0; dr < nDr; ++dr)
- {
-// const size_t drIdx = drUp[dr];
-// if(drIdx >= nDrRes)
- if(dr >= nDrRes)
- {
- break;
- }
-
- // Re-simulate for residual if required.
- if(timeFactor != 1 || nCh != nChRes)
- {
- fill(&(residual[thread][0][0][0]),
- &(residual[thread][0][0][0]) + nBl * nChRes * nCr,
- dcomplex());
-
- const_cursor<double> cr_uvw = casa_const_cursor(target[i
- * timeFactor + j].getUVW());
- splitUVW(nSt, nBl, cr_baseline, cr_uvw, cr_split);
-
-// rotateUVW(itsPhaseRef, itsPatchList[drIdx]->position(),
-// nSt, cr_split);
- rotateUVW(itsPhaseRef, itsPatchList[dr]->position(),
- nSt, cr_split);
-
-// simulate(itsPatchList[drIdx]->position(),
-// itsPatchList[drIdx], nSt, nBl, nChRes, cr_baseline,
-// cr_freqRes, cr_split, cr_model_res);
- simulate(itsPatchList[dr]->position(), itsPatchList[dr],
- nSt, nBl, nChRes, cr_baseline, cr_freqRes, cr_split,
- cr_model_res);
- }
- else
- {
- copy(&(buffer[thread][dr][0][0][0]),
- &(buffer[thread][dr][0][0][0]) + nBl * nChRes * nCr,
- &(residual[thread][0][0][0]));
- }
-
- // Apply solutions.
- size_t strides_coeff[2] = {1, 8};
-// const_cursor<double> cr_coeff(&(unknowns[thread][drIdx * nSt
-// * 8]), 2, strides_coeff);
- const_cursor<double> cr_coeff(&(unknowns[thread][dr * nSt
- * 8]), 2, strides_coeff);
- apply(nBl, nChRes, cr_baseline, cr_coeff, cr_model_res);
-
- // Subtract.
- cursor<fcomplex> cr_residual =
- casa_cursor(target[i * timeFactor + j].getData());
-
- const IPosition tmp_strides_mix_res =
- itsFactorsSubtr[i * timeFactor + j].steps();
-
- // The target direction is always last, and therefore it has
- // index itsNDir - 1. The directions to subtract are always
- // first, and therefore have indices [0, nDrRes).
-// const_cursor<dcomplex> cr_mix(&(itsFactorsSubtr[i
-// * timeFactor + j](IPosition(5, itsNDir - 1, drIdx, 0, 0,
-// 0))), 3, tmp_strides_mix_res.storage() + 2);
- const_cursor<dcomplex> cr_mix(&(itsFactorsSubtr[i
- * timeFactor + j](IPosition(5, itsNDir - 1, dr, 0, 0,
- 0))), 3, tmp_strides_mix_res.storage() + 2);
- subtract(nBl, nChRes, cr_baseline, cr_residual,
- cr_model_res, cr_mix);
- }
- }
+ // Estimate Jones matrices.
+ //
+ // A Jones matrix will be estimated for each pair of station and
+ // direction.
+ //
+ // A single (overdetermined) non-linear set of equations for all
+ // stations and directions is solved iteratively. The influence of
+ // each direction on each other direction is given by the mixing
+ // matrix.
+ const_cursor<bool> cr_flag =
+ casa_const_cursor(itsAvgResults[0]->get()[ts].getFlags());
+ const_cursor<float> cr_weight =
+ casa_const_cursor(itsAvgResults[0]->get()[ts].getWeights());
+ const_cursor<dcomplex> cr_mix = casa_const_cursor(itsFactors[ts]);
+
+ vector<const_cursor<fcomplex> > cr_data(nDr);
+ vector<const_cursor<dcomplex> > cr_model(nDr);
+ for(size_t dr = 0; dr < nDr; ++dr)
+ {
+ cr_data[dr] =
+ casa_const_cursor(itsAvgResults[dr]->get()[ts].getData());
+ cr_model[dr] =
+ const_cursor<dcomplex>(&(storage.model[dr * nSamples]), 3,
+ stride_model);
+ }
- // Copy solutions to global solution array.
-// for(size_t dr = 0; dr < nDrUp; ++dr)
- for(size_t dr = 0; dr < nDr; ++dr)
- {
-// size_t idx = drUp[dr];
-// copy(&(unknowns[thread][idx * nSt * 8]), &(unknowns[thread][idx
-// * nSt * 8]) + nSt * 8, &(itsUnknowns(IPosition(4, 0, 0, idx,
-// itsTimeCount + i))));
-// copy(&(errors[thread][idx * nSt * 8]), &(errors[thread][idx
-// * nSt * 8]) + nSt * 8, &(itsErrors(IPosition(4, 0, 0, idx,
-// itsTimeCount + i))));
- copy(&(unknowns[thread][dr * nSt * 8]), &(unknowns[thread][dr
- * nSt * 8]) + nSt * 8, &(itsUnknowns(IPosition(4, 0, 0, dr,
- itsTimeCount + i))));
- }
+ bool converged = estimate(nDr, nSt, nBl, nCh, cr_baseline, cr_data,
+ cr_model, cr_flag, cr_weight, cr_mix, &(storage.unknowns[0]));
+ if(converged)
+ {
+ ++storage.count_converged;
}
- // Store last known solutions.
- for(size_t dr = 0; dr < nDr; ++dr)
+ // Compute the residual.
+ //
+ // All the so-called "subtract sources" are subtracted from the
+ // observed data. The previously estimated Jones matrices, as well as
+ // the appropriate mixing weight are applied before subtraction.
+ //
+ // Note that the resolution of the residual can differ from the
+ // resolution at which the Jones matrices were estimated.
+ for(size_t ts_subtr = multiplier * ts, ts_subtr_end = min(ts_subtr
+ + multiplier, nTimeSubtr); ts_subtr != ts_subtr_end; ++ts_subtr)
{
-// int idx = last[dr][0];
-// for(size_t i = 1; i < nThread; ++i)
-// {
-// idx = max(idx, last[dr][i]);
-// }
- int idx = (nTime > 0 ? (nTime - 1) : -1);
-
- if(idx >= 0)
+ for(size_t dr = 0; dr < nDrSubtr; ++dr)
+ {
+ // Re-use simulation used for estimating Jones matrices if possible.
+ cursor<dcomplex> cr_model_subtr(&(storage.model[dr * nSamples]),
+ 3, stride_model);
+
+ // Re-simulate if required.
+ if(multiplier != 1 || nCh != nChSubtr)
{
- copy(&(itsUnknowns(IPosition(4, 0, 0, dr, itsTimeCount + idx))),
- &(itsUnknowns(IPosition(4, 0, 0, dr, itsTimeCount + idx)))
- + nSt * 8, &(itsLastKnowns(IPosition(3, 0, 0, dr))));
+ const_cursor<double> cr_uvw =
+ casa_const_cursor(itsAvgResultSubtr->get()[ts_subtr].getUVW());
+ splitUVW(nSt, nBl, cr_baseline, cr_uvw, cr_uvw_split);
+
+ // Rotate the UVW coordinates for the target direction to the
+ // direction of source to subtract. This is required because at
+ // the resolution of the residual the UVW coordinates for
+ // directions other than the target are unavailable (unless the
+ // resolution of the residual is equal to the resolution at which
+ // the Jones matrices were estimated, of course).
+ rotateUVW(itsPhaseRef, itsPatchList[dr]->position(), nSt,
+ cr_uvw_split);
+
+ // Zero the visibility buffer.
+ fill(storage.model_subtr.begin(), storage.model_subtr.end(),
+ dcomplex());
+
+ // Simulate visibilities at the resolution of the residual.
+ size_t stride_model_subtr[3] = {1, nCr, nCr * nChSubtr};
+ cr_model_subtr = cursor<dcomplex>(&(storage.model_subtr[0]), 3,
+ stride_model_subtr);
+ simulate(itsPatchList[dr]->position(), itsPatchList[dr], nSt, nBl,
+ nChSubtr, cr_baseline, cr_freqSubtr, cr_uvw_split,
+ cr_model_subtr);
}
- }
- // Update convergence count.
- for (size_t i=0; i<nThread; ++i) {
- itsNConverged += converged[i];
+ // Apply Jones matrices.
+ size_t stride_unknowns[2] = {1, 8};
+ const_cursor<double> cr_unknowns(&(storage.unknowns[dr * nSt * 8]),
+ 2, stride_unknowns);
+ apply(nBl, nChSubtr, cr_baseline, cr_unknowns, cr_model_subtr);
+
+ // Subtract the source contribution from the data.
+ cursor<fcomplex> cr_residual =
+ casa_cursor(itsAvgResultSubtr->get()[ts_subtr].getData());
+
+ // Construct a cursor to iterate over a slice of the mixing matrix
+ // at the resolution of the residual. The "to" and "from" direction
+ // are fixed. Since the full mixing matrix is 5-D, the slice is
+ // therefore 3-D. Each individual value in the slice quantifies the
+ // influence of the source to subtract on the target direction for
+ // a particular correlation, channel, and baseline.
+ //
+ // The target direction is the direction with the highest index by
+ // convention, i.e. index itsNDir - 1. The directions to subtract
+ // have the lowest indices by convention, i.e. indices
+ // [0, nDrSubtr).
+ const IPosition &stride_mix_subtr =
+ itsFactorsSubtr[ts_subtr].steps();
+ size_t stride_mix_subtr_slice[3] = {stride_mix_subtr[2],
+ stride_mix_subtr[3], stride_mix_subtr[4]};
+ ASSERT(stride_mix_subtr_slice[0] == itsNDir * itsNDir
+ && stride_mix_subtr_slice[1] == itsNDir * itsNDir * nCr
+ && stride_mix_subtr_slice[2] == itsNDir * itsNDir * nCr * nChSubtr);
+
+ IPosition offset(5, itsNDir - 1, dr, 0, 0, 0);
+ const_cursor<dcomplex> cr_mix_subtr =
+ const_cursor<dcomplex>(&(itsFactorsSubtr[ts_subtr](offset)), 3,
+ stride_mix_subtr_slice);
+
+ // Subtract the source.
+ subtract(nBl, nChSubtr, cr_baseline, cr_residual, cr_model_subtr,
+ cr_mix_subtr);
+ }
}
- itsTimerSolve.stop();
+ // Copy solutions to global solution array.
+ copy(storage.unknowns.begin(), storage.unknowns.end(),
+ &(itsUnknowns[(itsTimeIndex + ts) * nDr * nSt * 8]));
}
- // Let the next step process the data.
- itsTimer.stop();
- for (uint i=0; i<itsNTimeOutSubtr; ++i) {
- getNextStep()->process (itsAvgResultSubtr->get()[i]);
- itsAvgResultSubtr->get()[i].clear();
+ // Store last known solutions.
+ if(nTime > 0)
+ {
+ copy(&(itsUnknowns[(itsTimeIndex + nTime - 1) * nDr * nSt * 8]),
+ &(itsUnknowns[(itsTimeIndex + nTime) * nDr * nSt * 8]),
+ itsLastKnowns.begin());
}
- // Clear the vector in the MultiStep.
- itsAvgResultSubtr->clear();
- itsTimer.start();
- }
-
- string Demixer::toString (double value) const
- {
- ostringstream os;
- os << setprecision(16) << value;
- return os.str();
- }
- double Demixer::getAngle (const String& value) const
- {
- double angle;
- Quantity q;
- ASSERTSTR (Quantity::read (q, value),
- "Demixer: " + value + " is not a proper angle");
- if (q.getUnit().empty()) {
- angle = q.getValue() / 180. * C::pi;
- } else {
- ASSERTSTR (q.getFullUnit().getValue() == UnitVal::ANGLE,
- "Demixer: " + value + " is not a proper angle");
- angle = q.getValue("rad");
+ // Update convergence count.
+ for(size_t i = 0; i < nThread; ++i)
+ {
+ itsNConverged += threadStorage[i].count_converged;
}
- return angle;
}
void Demixer::dumpSolutions()
@@ -1125,10 +1012,8 @@ namespace LOFAR {
// Write solutions.
for(size_t ts = 0; ts < itsNTimeDemix; ++ts) {
- double *unknowns = &(itsUnknowns(IPosition(4, 0, 0, 0, ts)));
-// double *errors = &(itsErrors(IPosition(4, 0, 0, 0, ts)));
+ double *unknowns = &(itsUnknowns[ts * itsNModel * itsNStation * 8]);
for(size_t i = 0; i < parms.size(); ++i) {
-// parms[i].setCoeff(BBS::Location(0, ts), unknowns + i, 1, errors + i);
parms[i].setCoeff(BBS::Location(0, ts), unknowns + i, 1);
}
}
@@ -1137,6 +1022,31 @@ namespace LOFAR {
parmCache.flush();
}
+ namespace
+ {
+ string toString (double value)
+ {
+ ostringstream os;
+ os << setprecision(16) << value;
+ return os.str();
+ }
+
+ double getAngle (const String& value)
+ {
+ double angle;
+ Quantity q;
+ ASSERTSTR (Quantity::read (q, value),
+ "Demixer: " + value + " is not a proper angle");
+ if (q.getUnit().empty()) {
+ angle = q.getValue() / 180. * C::pi;
+ } else {
+ ASSERTSTR (q.getFullUnit().getValue() == UnitVal::ANGLE,
+ "Demixer: " + value + " is not a proper angle");
+ angle = q.getValue("rad");
+ }
+ return angle;
+ }
+ } //# end unnamed namespace
} //# end namespace DPPP
} //# end namespace LOFAR
diff --git a/CEP/DP3/DPPP/src/EstimateMixed.cc b/CEP/DP3/DPPP/src/EstimateMixed.cc
index 091d4773077..2cdbffa5af3 100644
--- a/CEP/DP3/DPPP/src/EstimateMixed.cc
+++ b/CEP/DP3/DPPP/src/EstimateMixed.cc
@@ -27,18 +27,50 @@
#include <DPPP/EstimateMixed.h>
#include <Common/LofarLogger.h>
#include <scimath/Fitting/LSQFit.h>
-#include <boost/multi_array.hpp>
namespace LOFAR
{
namespace DPPP
{
+namespace
+{
+// Compute a map that contains the index of the unknowns related to the
+// specified baseline in the list of all unknowns.
+void makeIndex(size_t nDirection, size_t nStation, const Baseline &baseline,
+ unsigned int *index)
+{
+ const size_t nCorrelation = 4;
+ for(size_t cr = 0; cr < nCorrelation; ++cr)
+ {
+ size_t idx0 = baseline.first * 8 + (cr / 2) * 4;
+ size_t idx1 = baseline.second * 8 + (cr % 2) * 4;
+
+ for(size_t dr = 0; dr < nDirection; ++dr)
+ {
+ *index++ = idx0;
+ *index++ = idx0 + 1;
+ *index++ = idx0 + 2;
+ *index++ = idx0 + 3;
+
+ *index++ = idx1;
+ *index++ = idx1 + 1;
+ *index++ = idx1 + 2;
+ *index++ = idx1 + 3;
+
+ idx0 += nStation * 8;
+ idx1 += nStation * 8;
+ }
+ }
+}
+} // Unnamed namespace.
+
+
bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
size_t nChannel, const_cursor<Baseline> baselines,
vector<const_cursor<fcomplex> > data, vector<const_cursor<dcomplex> > model,
const_cursor<bool> flag, const_cursor<float> weight,
- const_cursor<dcomplex> mix, double *unknowns, double *errors)
+ const_cursor<dcomplex> mix, double *unknowns)
{
ASSERT(data.size() == nDirection && model.size() == nDirection);
@@ -52,46 +84,13 @@ bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
// each visibility provides information about (no. of directions) x 2 x 2
// x 2 (scalar) unknowns = (no. of directions) x 8. For each of these
// unknowns the value of the partial derivative of the model with respect
- // to the unknow has to be computed.
- const unsigned int nPartial = nDirection * 8;
-
- // Construct partial derivative index template for each correlation.
- boost::multi_array<unsigned int, 2> dIndexTemplate(boost::extents[4]
- [nPartial]);
- for(size_t cr = 0; cr < 4; ++cr)
- {
- size_t idx0 = (cr / 2) * 4;
- size_t idx1 = (cr & 1) * 4;
-
- for(size_t dr = 0; dr < nDirection; ++dr)
- {
- dIndexTemplate[cr][dr * 8 + 0] = idx0 + 0;
- dIndexTemplate[cr][dr * 8 + 1] = idx0 + 1;
- dIndexTemplate[cr][dr * 8 + 2] = idx0 + 2;
- dIndexTemplate[cr][dr * 8 + 3] = idx0 + 3;
- dIndexTemplate[cr][dr * 8 + 4] = idx1 + 0;
- dIndexTemplate[cr][dr * 8 + 5] = idx1 + 1;
- dIndexTemplate[cr][dr * 8 + 6] = idx1 + 2;
- dIndexTemplate[cr][dr * 8 + 7] = idx1 + 3;
- idx0 += nStation * 8;
- idx1 += nStation * 8;
- }
- }
+ // to the unknown has to be computed.
+ const size_t nPartial = nDirection * 8;
+ vector<unsigned int> dIndex(4 * nPartial);
// Allocate space for intermediate results.
- boost::multi_array<dcomplex, 2> M(boost::extents[nDirection][4]);
- boost::multi_array<dcomplex, 3> dM(boost::extents[nDirection][4][4]);
- boost::multi_array<double, 1> dR(boost::extents[nPartial]);
- boost::multi_array<double, 1> dI(boost::extents[nPartial]);
- boost::multi_array<unsigned int, 2> dIndex(boost::extents[4][nPartial]);
-
- dcomplex coherence;
- dcomplex Jp_00, Jp_01, Jp_10, Jp_11;
- dcomplex Jq_00, Jq_01, Jq_10, Jq_11;
- dcomplex Jp_00_s0, Jp_10_s0, Jp_00_s1, Jp_10_s1, Jp_01_s2, Jp_11_s2,
- Jp_01_s3, Jp_11_s3;
- dcomplex Jq_00_s0, Jq_10_s0, Jq_01_s1, Jq_11_s1, Jq_00_s2, Jq_10_s2,
- Jq_01_s3, Jq_11_s3;
+ vector<dcomplex> M(nDirection * 4), dM(nDirection * 16);
+ vector<double> dR(nPartial), dI(nPartial);
// Iterate until convergence.
size_t nIterations = 0;
@@ -104,103 +103,78 @@ bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
if(p != q)
{
- // Update partial derivative index for current baseline.
- const size_t offsetP = p * 8;
- const size_t offsetQ = q * 8;
- for(size_t cr = 0; cr < 4; ++cr)
- {
- for(size_t dr = 0; dr < nDirection; ++dr)
- {
- dIndex[cr][dr * 8 + 0] = dIndexTemplate[cr][dr * 8 + 0]
- + offsetP;
- dIndex[cr][dr * 8 + 1] = dIndexTemplate[cr][dr * 8 + 1]
- + offsetP;
- dIndex[cr][dr * 8 + 2] = dIndexTemplate[cr][dr * 8 + 2]
- + offsetP;
- dIndex[cr][dr * 8 + 3] = dIndexTemplate[cr][dr * 8 + 3]
- + offsetP;
-
- dIndex[cr][dr * 8 + 4] = dIndexTemplate[cr][dr * 8 + 4]
- + offsetQ;
- dIndex[cr][dr * 8 + 5] = dIndexTemplate[cr][dr * 8 + 5]
- + offsetQ;
- dIndex[cr][dr * 8 + 6] = dIndexTemplate[cr][dr * 8 + 6]
- + offsetQ;
- dIndex[cr][dr * 8 + 7] = dIndexTemplate[cr][dr * 8 + 7]
- + offsetQ;
- }
- }
+ // Create partial derivative index for current baseline.
+ makeIndex(nDirection, nStation, *baselines, &(dIndex[0]));
for(size_t ch = 0; ch < nChannel; ++ch)
{
- dcomplex coherence;
for(size_t dr = 0; dr < nDirection; ++dr)
{
+ // Jones matrix for station P.
const double *Jp =
&(unknowns[dr * nStation * 8 + p * 8]);
- Jp_00 = dcomplex(Jp[0], Jp[1]);
- Jp_01 = dcomplex(Jp[2], Jp[3]);
- Jp_10 = dcomplex(Jp[4], Jp[5]);
- Jp_11 = dcomplex(Jp[6], Jp[7]);
+ const dcomplex Jp_00(Jp[0], Jp[1]);
+ const dcomplex Jp_01(Jp[2], Jp[3]);
+ const dcomplex Jp_10(Jp[4], Jp[5]);
+ const dcomplex Jp_11(Jp[6], Jp[7]);
+ // Jones matrix for station Q, conjugated.
const double *Jq =
&(unknowns[dr * nStation * 8 + q * 8]);
- Jq_00 = dcomplex(Jq[0], -Jq[1]);
- Jq_01 = dcomplex(Jq[2], -Jq[3]);
- Jq_10 = dcomplex(Jq[4], -Jq[5]);
- Jq_11 = dcomplex(Jq[6], -Jq[7]);
-
- coherence = (model[dr][0]);
- Jp_00_s0 = Jp_00 * coherence;
- Jp_10_s0 = Jp_10 * coherence;
- Jq_00_s0 = Jq_00 * coherence;
- Jq_10_s0 = Jq_10 * coherence;
-
- coherence = (model[dr][1]);
- Jp_00_s1 = Jp_00 * coherence;
- Jp_10_s1 = Jp_10 * coherence;
- Jq_01_s1 = Jq_01 * coherence;
- Jq_11_s1 = Jq_11 * coherence;
-
- coherence = (model[dr][2]);
- Jp_01_s2 = Jp_01 * coherence;
- Jp_11_s2 = Jp_11 * coherence;
- Jq_00_s2 = Jq_00 * coherence;
- Jq_10_s2 = Jq_10 * coherence;
-
- coherence = (model[dr][3]);
- Jp_01_s3 = Jp_01 * coherence;
- Jp_11_s3 = Jp_11 * coherence;
- Jq_01_s3 = Jq_01 * coherence;
- Jq_11_s3 = Jq_11 * coherence;
-
- M[dr][0] = Jp_00 * (Jq_00_s0 + Jq_01_s1)
- + Jp_01 * (Jq_00_s2 + Jq_01_s3);
- dM[dr][0][0] = Jq_00_s0 + Jq_01_s1;
- dM[dr][0][1] = Jq_00_s2 + Jq_01_s3;
- dM[dr][0][2] = Jp_00_s0 + Jp_01_s2;
- dM[dr][0][3] = Jp_00_s1 + Jp_01_s3;
-
- M[dr][1] = Jp_00 * (Jq_10_s0 + Jq_11_s1)
- + Jp_01 * (Jq_10_s2 + Jq_11_s3);
- dM[dr][1][0] = Jq_10_s0 + Jq_11_s1;
- dM[dr][1][1] = Jq_10_s2 + Jq_11_s3;
- dM[dr][1][2] = dM[dr][0][2];
- dM[dr][1][3] = dM[dr][0][3];
-
- M[dr][2] = Jp_10 * (Jq_00_s0 + Jq_01_s1)
- + Jp_11 * (Jq_00_s2 + Jq_01_s3);
- dM[dr][2][0] = dM[dr][0][0];
- dM[dr][2][1] = dM[dr][0][1];
- dM[dr][2][2] = Jp_10_s0 + Jp_11_s2;
- dM[dr][2][3] = Jp_10_s1 + Jp_11_s3;
-
- M[dr][3] = Jp_10 * (Jq_10_s0 + Jq_11_s1)
- + Jp_11 * (Jq_10_s2 + Jq_11_s3);
- dM[dr][3][0] = dM[dr][1][0];
- dM[dr][3][1] = dM[dr][1][1];
- dM[dr][3][2] = dM[dr][2][2];
- dM[dr][3][3] = dM[dr][2][3];
+ const dcomplex Jq_00(Jq[0], -Jq[1]);
+ const dcomplex Jq_01(Jq[2], -Jq[3]);
+ const dcomplex Jq_10(Jq[4], -Jq[5]);
+ const dcomplex Jq_11(Jq[6], -Jq[7]);
+
+ // Fetch model visibilities for the current direction.
+ const dcomplex xx = model[dr][0];
+ const dcomplex xy = model[dr][1];
+ const dcomplex yx = model[dr][2];
+ const dcomplex yy = model[dr][3];
+
+ // Precompute terms involving conj(Jq) and the model
+ // visibilities.
+ const dcomplex Jq_00xx_01xy = Jq_00 * xx + Jq_01 * xy;
+ const dcomplex Jq_00yx_01yy = Jq_00 * yx + Jq_01 * yy;
+ const dcomplex Jq_10xx_11xy = Jq_10 * xx + Jq_11 * xy;
+ const dcomplex Jq_10yx_11yy = Jq_10 * yx + Jq_11 * yy;
+
+ // Precompute (Jp x conj(Jq)) * vec(data), where 'x'
+ // denotes the Kronecker product. This is the model
+ // visibility for the current direction, with the
+ // current Jones matrix estimates applied. This is
+ // stored in M.
+ // Also, precompute the partial derivatives of M with
+ // respect to all 16 parameters (i.e. 2 Jones matrices
+ // Jp and Jq, 4 complex scalars per Jones matrix, 2 real
+ // scalars per complex scalar, 2 * 4 * 2 = 16). These
+ // partial derivatives are stored in dM.
+ M[dr * 4] = Jp_00 * Jq_00xx_01xy + Jp_01 * Jq_00yx_01yy;
+ dM[dr * 16] = Jq_00xx_01xy;
+ dM[dr * 16 + 1] = Jq_00yx_01yy;
+ dM[dr * 16 + 2] = Jp_00 * xx + Jp_01 * yx;
+ dM[dr * 16 + 3] = Jp_00 * xy + Jp_01 * yy;
+
+ M[dr * 4 + 1] = Jp_00 * Jq_10xx_11xy + Jp_01
+ * Jq_10yx_11yy;
+ dM[dr * 16 + 4] = Jq_10xx_11xy;
+ dM[dr * 16 + 5] = Jq_10yx_11yy;
+ dM[dr * 16 + 6] = dM[dr * 16 + 2];
+ dM[dr * 16 + 7] = dM[dr * 16 + 3];
+
+ M[dr * 4 + 2] = Jp_10 * Jq_00xx_01xy + Jp_11
+ * Jq_00yx_01yy;
+ dM[dr * 16 + 8] = dM[dr * 16];
+ dM[dr * 16 + 9] = dM[dr * 16 + 1];
+ dM[dr * 16 + 10] = Jp_10 * xx + Jp_11 * yx;
+ dM[dr * 16 + 11] = Jp_10 * xy + Jp_11 * yy;
+
+ M[dr * 4 + 3] = Jp_10 * Jq_10xx_11xy + Jp_11
+ * Jq_10yx_11yy;
+ dM[dr * 16 + 12] = dM[dr * 16 + 4];
+ dM[dr * 16 + 13] = dM[dr * 16 + 5];
+ dM[dr * 16 + 14] = dM[dr * 16 + 10];
+ dM[dr * 16 + 15] = dM[dr * 16 + 11];
}
for(size_t cr = 0; cr < 4; ++cr)
@@ -209,39 +183,44 @@ bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
{
for(size_t tg = 0; tg < nDirection; ++tg)
{
- dcomplex model = 0.0, partial;
+ dcomplex visibility(0.0, 0.0);
for(size_t dr = 0; dr < nDirection; ++dr)
{
- // Look-up mixing term.
- const dcomplex term = *mix;
-
- // Update model visibility.
- model += term * M[dr][cr];
-
- // Compute partial derivatives.
- partial = term * dM[dr][cr][0];
- dR[dr * 8] = real(partial);
- dI[dr * 8] = imag(partial);
- dR[dr * 8 + 1] = -imag(partial);
- dI[dr * 8 + 1] = real(partial);
-
- partial = term * dM[dr][cr][1];
- dR[dr * 8 + 2] = real(partial);
- dI[dr * 8 + 2] = imag(partial);
- dR[dr * 8 + 3] = -imag(partial);
- dI[dr * 8 + 3] = real(partial);
-
- partial = term * dM[dr][cr][2];
- dR[dr * 8 + 4] = real(partial);
- dI[dr * 8 + 4] = imag(partial);
- dR[dr * 8 + 5] = imag(partial);
- dI[dr * 8 + 5] = -real(partial);
-
- partial = term * dM[dr][cr][3];
- dR[dr * 8 + 6] = real(partial);
- dI[dr * 8 + 6] = imag(partial);
- dR[dr * 8 + 7] = imag(partial);
- dI[dr * 8 + 7] = -real(partial);
+ // Look-up mixing weight.
+ const dcomplex mix_weight = *mix;
+
+ // Weight model visibility.
+ visibility += mix_weight * M[dr * 4 + cr];
+
+ // Compute weighted partial derivatives.
+ dcomplex derivative(0.0, 0.0);
+ derivative =
+ mix_weight * dM[dr * 16 + cr * 4];
+ dR[dr * 8] = real(derivative);
+ dI[dr * 8] = imag(derivative);
+ dR[dr * 8 + 1] = -imag(derivative);
+ dI[dr * 8 + 1] = real(derivative);
+
+ derivative =
+ mix_weight * dM[dr * 16 + cr * 4 + 1];
+ dR[dr * 8 + 2] = real(derivative);
+ dI[dr * 8 + 2] = imag(derivative);
+ dR[dr * 8 + 3] = -imag(derivative);
+ dI[dr * 8 + 3] = real(derivative);
+
+ derivative =
+ mix_weight * dM[dr * 16 + cr * 4 + 2];
+ dR[dr * 8 + 4] = real(derivative);
+ dI[dr * 8 + 4] = imag(derivative);
+ dR[dr * 8 + 5] = imag(derivative);
+ dI[dr * 8 + 5] = -real(derivative);
+
+ derivative =
+ mix_weight * dM[dr * 16 + cr * 4 + 3];
+ dR[dr * 8 + 6] = real(derivative);
+ dI[dr * 8 + 6] = imag(derivative);
+ dR[dr * 8 + 7] = imag(derivative);
+ dI[dr * 8 + 7] = -real(derivative);
// Move to next source direction.
mix.forward(1);
@@ -249,14 +228,17 @@ bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
// Compute the residual.
dcomplex residual =
- static_cast<dcomplex>(data[tg][cr]) - model;
+ static_cast<dcomplex>(data[tg][cr])
+ - visibility;
// Update the normal equations.
- solver.makeNorm(nPartial, &(dIndex[cr][0]),
- &(dR[0]), static_cast<double>(weight[cr]),
+ solver.makeNorm(nPartial,
+ &(dIndex[cr * nPartial]), &(dR[0]),
+ static_cast<double>(weight[cr]),
real(residual));
- solver.makeNorm(nPartial, &(dIndex[cr][0]),
- &(dI[0]), static_cast<double>(weight[cr]),
+ solver.makeNorm(nPartial,
+ &(dIndex[cr * nPartial]), &(dI[0]),
+ static_cast<double>(weight[cr]),
imag(residual));
// Move to next target direction.
@@ -328,22 +310,6 @@ bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
++nIterations;
}
- // Get the estimated error for each unknown from the solver.
- if(errors)
- {
-// boost::multi_array<double, 1> errors_packed(boost::extents[nUnknowns * nUnknowns]);
-// // TODO: Somehow this only returns zero??
-// bool status = solver.getErrors(&(errors_packed[0]));
-
- vector<double> cov(nUnknowns * nUnknowns);
- bool status = solver.getCovariance(&(cov[0]));
- ASSERT(status);
- for(size_t i = 0; i < nUnknowns; ++i)
- {
- errors[i] = sqrt(abs(cov[i * nUnknowns + i])) * solver.getSD();
- }
- }
-
const bool converged = (solver.isReady() == casa::LSQFit::SOLINCREMENT
|| solver.isReady() == casa::LSQFit::DERIVLEVEL);
return converged;
diff --git a/CEP/DP3/DPPP/src/Patch.cc b/CEP/DP3/DPPP/src/Patch.cc
index 6280d598bec..d75984adebb 100644
--- a/CEP/DP3/DPPP/src/Patch.cc
+++ b/CEP/DP3/DPPP/src/Patch.cc
@@ -30,15 +30,6 @@ namespace LOFAR
namespace DPPP
{
-void Patch::accept(ModelComponentVisitor &visitor) const
-{
- for(const_iterator it = begin(), it_end = end(); it != it_end; ++it)
- {
- (*it)->accept(visitor);
- }
- visitor.visit(*this);
-}
-
Patch::const_iterator Patch::begin() const
{
return itsComponents.begin();
@@ -53,27 +44,25 @@ void Patch::computePosition()
{
itsPosition = Position();
- if(itsComponents.empty())
+ if(!itsComponents.empty())
{
- return;
- }
+ double x = 0.0, y = 0.0, z = 0.0;
+ for(const_iterator it = begin(), it_end = end(); it != it_end; ++it)
+ {
+ const Position &position = (*it)->position();
+ double cosDec = cos(position[1]);
+ x += cos(position[0]) * cosDec;
+ y += sin(position[0]) * cosDec;
+ z += sin(position[1]);
+ }
- double x = 0.0, y = 0.0, z = 0.0;
- for(const_iterator it = begin(), it_end = end(); it != it_end; ++it)
- {
- const Position &position = (*it)->position();
- double cosDec = cos(position[1]);
- x += cos(position[0]) * cosDec;
- y += sin(position[0]) * cosDec;
- z += sin(position[1]);
- }
+ x /= itsComponents.size();
+ y /= itsComponents.size();
+ z /= itsComponents.size();
- x /= itsComponents.size();
- y /= itsComponents.size();
- z /= itsComponents.size();
-
- itsPosition[0] = atan2(y, x);
- itsPosition[1] = asin(z);
+ itsPosition[0] = atan2(y, x);
+ itsPosition[1] = asin(z);
+ }
}
} //# namespace DPPP
diff --git a/CEP/DP3/DPPP/src/PhaseShift.cc b/CEP/DP3/DPPP/src/PhaseShift.cc
index 70610f65e91..c2475222056 100644
--- a/CEP/DP3/DPPP/src/PhaseShift.cc
+++ b/CEP/DP3/DPPP/src/PhaseShift.cc
@@ -137,30 +137,30 @@ namespace LOFAR {
//# to process.
#pragma omp parallel for
for (int i=0; i<nbl; ++i) {
- Complex* __restrict__ data = newBuf.getData().data() + i*nchan*ncorr;
- double* __restrict__ uvw = newBuf.getUVW().data() + i*3;
- DComplex* __restrict__ phasors = itsPhasors.data() + i*nchan;
- double u = uvw[0]*mat1[0] + uvw[1]*mat1[3] + uvw[2]*mat1[6];
- double v = uvw[0]*mat1[1] + uvw[1]*mat1[4] + uvw[2]*mat1[7];
- double w = uvw[0]*mat1[2] + uvw[1]*mat1[5] + uvw[2]*mat1[8];
- double phase = itsXYZ[0]*uvw[0] + itsXYZ[1]*uvw[1] + itsXYZ[2]*uvw[2];
- for (int j=0; j<nchan; ++j) {
- // Shift the phase of the data of this baseline.
- // Converting the phase term to wavelengths (and applying 2*pi)
- // u_wvl = u_m / wvl = u_m * freq / c
- // has been done once in the beginning (in updateInfo).
- double phasewvl = phase * itsFreqC[j];
- DComplex phasor(cos(phasewvl), sin(phasewvl));
- *phasors++ = phasor;
- for (int k=0; k<ncorr; ++k) {
- *data = DComplex(*data) * phasor;
- data++;
- }
- }
- uvw[0] = u;
- uvw[1] = v;
- uvw[2] = w;
- uvw += 3;
+ Complex* __restrict__ data = newBuf.getData().data() + i*nchan*ncorr;
+ double* __restrict__ uvw = newBuf.getUVW().data() + i*3;
+ DComplex* __restrict__ phasors = itsPhasors.data() + i*nchan;
+ double u = uvw[0]*mat1[0] + uvw[1]*mat1[3] + uvw[2]*mat1[6];
+ double v = uvw[0]*mat1[1] + uvw[1]*mat1[4] + uvw[2]*mat1[7];
+ double w = uvw[0]*mat1[2] + uvw[1]*mat1[5] + uvw[2]*mat1[8];
+ double phase = itsXYZ[0]*uvw[0] + itsXYZ[1]*uvw[1] + itsXYZ[2]*uvw[2];
+ for (int j=0; j<nchan; ++j) {
+ // Shift the phase of the data of this baseline.
+ // Converting the phase term to wavelengths (and applying 2*pi)
+ // u_wvl = u_m / wvl = u_m * freq / c
+ // has been done once in the beginning (in updateInfo).
+ double phasewvl = phase * itsFreqC[j];
+ DComplex phasor(cos(phasewvl), sin(phasewvl));
+ *phasors++ = phasor;
+ for (int k=0; k<ncorr; ++k) {
+ *data = DComplex(*data) * phasor;
+ data++;
+ }
+ }
+ uvw[0] = u;
+ uvw[1] = v;
+ uvw[2] = w;
+ uvw += 3;
} //# end omp parallel for
itsTimer.stop();
getNextStep()->process (newBuf);
diff --git a/CEP/DP3/DPPP/src/PointSource.cc b/CEP/DP3/DPPP/src/PointSource.cc
index 5e6468257fa..6cfb47dfdda 100644
--- a/CEP/DP3/DPPP/src/PointSource.cc
+++ b/CEP/DP3/DPPP/src/PointSource.cc
@@ -37,7 +37,8 @@ PointSource::PointSource(const Position &position)
itsRefFreq(0.0),
itsPolarizedFraction(0.0),
itsPolarizationAngle(0.0),
- itsRotationMeasure(0.0)
+ itsRotationMeasure(0.0),
+ itsHasRotationMeasure(false)
{
}
@@ -47,7 +48,8 @@ PointSource::PointSource(const Position &position, const Stokes &stokes)
itsRefFreq(0.0),
itsPolarizedFraction(0.0),
itsPolarizationAngle(0.0),
- itsRotationMeasure(0.0)
+ itsRotationMeasure(0.0),
+ itsHasRotationMeasure(false)
{
}
@@ -61,19 +63,12 @@ void PointSource::setStokes(const Stokes &stokes)
itsStokes = stokes;
}
-void PointSource::setPolarizedFraction(double fraction)
+void PointSource::setRotationMeasure(double fraction, double angle, double rm)
{
itsPolarizedFraction = fraction;
-}
-
-void PointSource::setPolarizationAngle(double angle)
-{
itsPolarizationAngle = angle;
-}
-
-void PointSource::setRotationMeasure(double rm)
-{
itsRotationMeasure = rm;
+ itsHasRotationMeasure = true;
}
Stokes PointSource::stokes(double freq) const
@@ -124,12 +119,12 @@ void PointSource::accept(ModelComponentVisitor &visitor) const
bool PointSource::hasSpectralIndex() const
{
- return itsSpectralIndex.size() > 0;
+ return !itsSpectralIndex.empty();
}
bool PointSource::hasRotationMeasure() const
{
- return itsRotationMeasure > 0.0;
+ return itsHasRotationMeasure;
}
} //# namespace DPPP
diff --git a/CEP/DP3/DPPP/src/Simulate.cc b/CEP/DP3/DPPP/src/Simulate.cc
index 7b05559b346..200515705a5 100644
--- a/CEP/DP3/DPPP/src/Simulate.cc
+++ b/CEP/DP3/DPPP/src/Simulate.cc
@@ -22,12 +22,8 @@
#include <lofar_config.h>
#include <DPPP/Simulate.h>
-#include <DPPP/GaussianSource.h>
-#include <DPPP/ModelComponentVisitor.h>
-#include <DPPP/Patch.h>
-#include <DPPP/PointSource.h>
+#include <DPPP/Simulator.h>
#include <Common/LofarLogger.h>
-#include <casa/BasicSL/Constants.h>
// Only required for rotateUVW().
#include <DPPP/PhaseShift.h>
@@ -39,78 +35,165 @@ namespace LOFAR
namespace DPPP
{
-namespace
-{
-// Compute LMN coordinates of \p position relative to \p reference.
-//
-// \param[in] reference
-// Reference position on the celestial sphere.
-// \param[in] position
-// Position of interest on the celestial sphere.
-// \param[in] lmn
-// Pointer to a buffer of (at least) length three into which the computed LMN
-// coordinates will be written.
-void radec2lmn(const Position &reference, const Position &position,
- cursor<double> lmn);
-
-void phases(size_t nStation, size_t nChannel, const_cursor<double> lmn,
- const_cursor<double> uvw, const_cursor<double> freq,
- cursor<dcomplex> shift);
-
-void spectrum(const PointSource &source, size_t nChannel,
- const_cursor<double> freq, cursor<dcomplex> spectrum);
-} // Unnamed namespace.
-
void splitUVW(size_t nStation, size_t nBaseline,
const_cursor<Baseline> baselines, const_cursor<double> uvw,
cursor<double> split)
{
- cursor<double> known(split);
- vector<bool> flag(nStation, false);
-
- split[0] = 0.0;
- split[1] = 0.0;
- split[2] = 0.0;
- flag[0] = true;
-
+ // If the number of stations is zero, then the number of baselines should be
+ // zero as well because no valid baselines exist in this case.
+ ASSERT(nStation > 0 || nBaseline == 0);
+
+ // Flags to keep track of the stations for which the UVW coordinates are
+ // known.
+ vector<bool> flag(nStation, true);
+
+ // Find reachable stations, i.e. stations that participate in at least one
+ // (cross) baseline.
+ size_t nReachable = 0;
+ const_cursor<Baseline> tmp(baselines);
for(size_t i = 0; i < nBaseline; ++i)
{
- const size_t p = baselines->first;
- const size_t q = baselines->second;
- if(p != q && flag[p] != flag[q])
+ const size_t p = tmp->first;
+ const size_t q = tmp->second;
+
+ if(p != q)
{
if(flag[p])
{
- known.forward(1, p);
- split.forward(1, q);
- split[0] = uvw[0] + known[0];
- split[1] = uvw[1] + known[1];
- split[2] = uvw[2] + known[2];
- split.backward(1, q);
- known.backward(1, p);
- flag[q] = true;
+ flag[p] = false;
+ ++nReachable;
}
- else
+
+ if(flag[q])
{
- known.forward(1, q);
- split.forward(1, p);
- split[0] = -uvw[0] + known[0];
- split[1] = -uvw[1] + known[1];
- split[2] = -uvw[2] + known[2];
- split.backward(1, p);
- known.backward(1, q);
- flag[p] = true;
+ flag[q] = false;
+ ++nReachable;
}
}
+ if(nReachable == nStation)
+ {
+ break;
+ }
+
// Move to next baseline.
- uvw.forward(1);
- ++baselines;
- } // Baselines.
+ ++tmp;
+ }
- ASSERTSTR(static_cast<size_t>(std::count(flag.begin(), flag.end(), true))
- == flag.size(), "Unable to split baseline UVW coordinates into station"
- " UVW coordinates.");
+ // Zero the UVW coordinates of all unreachable stations and flag them as
+ // known.
+ if(nReachable < nStation)
+ {
+ for(size_t i = 0; i < nStation; ++i)
+ {
+ if(flag[i])
+ {
+ split.forward(1, i);
+ split[0] = 0.0;
+ split[1] = 0.0;
+ split[2] = 0.0;
+ split.backward(1, i);
+ }
+ }
+ }
+
+ // If no stations are reachable, nothing more can be done.
+ if(nReachable == 0)
+ {
+ return;
+ }
+
+ size_t ref = 0;
+ cursor<double> known(split);
+ while(true)
+ {
+ // Find first station for which UVW coordinates are unknown.
+ while(ref < nStation && flag[ref])
+ {
+ ++ref;
+ }
+
+ // UVW coordinates known for all stations, done.
+ if(ref == nStation)
+ {
+ break;
+ }
+
+ // Set UVW coordinates of the reference station to zero. Note that each
+ // isolated group of stations will have such a reference station. This
+ // is OK because the groups are isolated.
+ split.forward(1, ref);
+ split[0] = 0.0;
+ split[1] = 0.0;
+ split[2] = 0.0;
+ split.backward(1, ref);
+ flag[ref] = true;
+
+ // Single isolated station left, no use iterating over baselines.
+ if(ref == nStation - 1)
+ {
+ break;
+ }
+
+ bool done = false;
+ while(!done)
+ {
+ // Find UVW coordinates for other stations linked to the reference
+ // station via baselines.
+ size_t nFound = 0;
+ for(size_t i = 0; i < nBaseline; ++i)
+ {
+ const size_t p = baselines->first;
+ const size_t q = baselines->second;
+ if(p != q && flag[p] != flag[q])
+ {
+ if(flag[p])
+ {
+ known.forward(1, p);
+ split.forward(1, q);
+ split[0] = uvw[0] + known[0];
+ split[1] = uvw[1] + known[1];
+ split[2] = uvw[2] + known[2];
+ split.backward(1, q);
+ known.backward(1, p);
+ flag[q] = true;
+ }
+ else
+ {
+ known.forward(1, q);
+ split.forward(1, p);
+ split[0] = -uvw[0] + known[0];
+ split[1] = -uvw[1] + known[1];
+ split[2] = -uvw[2] + known[2];
+ split.backward(1, p);
+ known.backward(1, q);
+ flag[p] = true;
+ }
+
+ ++nFound;
+ }
+
+ // Move to next baseline.
+ uvw.forward(1);
+ ++baselines;
+ } // Baselines.
+
+ // Reset cursors.
+ uvw.backward(1, nBaseline);
+ baselines -= nBaseline;
+
+ // Depending on the baseline order it may be possible to derive the
+ // UVW coordinates of additional stations when UVW coordinates were
+ // found for at leat one station (i.e. nFound larger than 0).
+ // Another iteration is required in this case, unless UVW
+ // coordinates were found for all (remaining) stations (i.e. nFound
+ // equals nStation - ref - 1).
+ done = (nFound == 0 || nFound == nStation - ref - 1);
+ }
+ }
+
+ ASSERT(static_cast<size_t>(std::count(flag.begin(), flag.end(), true))
+ == flag.size());
}
void rotateUVW(const Position &from, const Position &to, size_t nUVW,
@@ -140,282 +223,18 @@ void rotateUVW(const Position &from, const Position &to, size_t nUVW,
} // Stations.
}
-class SimulateVisitor: public ModelComponentVisitor
-{
-public:
- SimulateVisitor(const Position &reference, size_t nStation,
- size_t nBaseline, size_t nChannel, const_cursor<Baseline> baselines,
- const_cursor<double> freq, const_cursor<double> uvw,
- cursor<dcomplex> buffer)
- : itsReference(reference),
- itsStationCount(nStation),
- itsBaselineCount(nBaseline),
- itsChannelCount(nChannel),
- itsBaselines(baselines),
- itsFreq(freq),
- itsUVW(uvw),
- itsBuffer(buffer),
- itsShiftBuffer(nStation * nChannel),
- itsSpectrumBuffer(nChannel * 4)
- {
- }
-
- void simulate(const ModelComponent::ConstPtr &component)
- {
- component->accept(*this);
- }
-
-private:
- virtual void visit(const PointSource &source)
- {
- // Compute LMN coordinates.
- double lmn[3];
- radec2lmn(itsReference, source.position(), cursor<double>(lmn));
-
- // Compute station phase shifts.
- phases(itsStationCount, itsChannelCount, const_cursor<double>(lmn),
- itsUVW, itsFreq, cursor<dcomplex>(&itsShiftBuffer[0]));
-
- // Compute source spectrum.
- spectrum(source, itsChannelCount, itsFreq,
- cursor<dcomplex>(&itsSpectrumBuffer[0]));
-
- // Compute visibilities.
- for(size_t bl = 0; bl < itsBaselineCount; ++bl)
- {
- const size_t p = itsBaselines->first;
- const size_t q = itsBaselines->second;
-
- if(p != q)
- {
- const dcomplex *shiftP = &(itsShiftBuffer[p * itsChannelCount]);
- const dcomplex *shiftQ = &(itsShiftBuffer[q * itsChannelCount]);
- const dcomplex *spectrum = &(itsSpectrumBuffer[0]);
- for(size_t ch = 0; ch < itsChannelCount; ++ch)
- {
- // Compute baseline phase shift.
- const dcomplex blShift = (*shiftQ) * conj(*shiftP);
- ++shiftP;
- ++shiftQ;
-
- // Compute visibilities.
- *itsBuffer += blShift * (*spectrum);
- ++itsBuffer;
- ++spectrum;
- *itsBuffer += blShift * (*spectrum);
- ++itsBuffer;
- ++spectrum;
- *itsBuffer += blShift * (*spectrum);
- ++itsBuffer;
- ++spectrum;
- *itsBuffer += blShift * (*spectrum);
- ++itsBuffer;
- ++spectrum;
-
- // Move to next channel.
- itsBuffer -= 4;
- itsBuffer.forward(1);
- } // Channels.
-
- itsBuffer.backward(1, itsChannelCount);
- }
-
- // Move to next baseline.
- itsBuffer.forward(2);
- ++itsBaselines;
- } // Baselines.
-
- itsBuffer.backward(2, itsBaselineCount);
- itsBaselines -= itsBaselineCount;
- }
-
- virtual void visit(const GaussianSource &source)
- {
- // Compute LMN coordinates.
- double lmn[3];
- radec2lmn(itsReference, source.position(), cursor<double>(lmn));
-
- // Compute station phase shifts.
- phases(itsStationCount, itsChannelCount, const_cursor<double>(lmn),
- itsUVW, itsFreq, cursor<dcomplex>(&itsShiftBuffer[0]));
-
- // Compute source spectrum.
- spectrum(source, itsChannelCount, itsFreq,
- cursor<dcomplex>(&itsSpectrumBuffer[0]));
-
- // Convert position angle from North over East to the angle used to
- // rotate the right-handed UV-plane.
- // TODO: Can probably optimize by changing the rotation matrix instead.
- const double phi = casa::C::pi_2 - source.positionAngle();
- const double cosPhi = cos(phi);
- const double sinPhi = sin(phi);
-
- // Take care of the conversion of axis lengths from FWHM in radians to
- // sigma.
- // TODO: Shouldn't the projection from the celestial sphere to the
- // UV-plane be taken into account here?
- const double fwhm2sigma = 1.0 / (2.0 * std::sqrt(2.0 * std::log(2.0)));
- const double uScale = source.majorAxis() * fwhm2sigma;
- const double vScale = source.minorAxis() * fwhm2sigma;
-
- for(size_t bl = 0; bl < itsBaselineCount; ++bl)
- {
- const size_t p = itsBaselines->first;
- const size_t q = itsBaselines->second;
-
- if(p != q)
- {
- itsUVW.forward(1, q);
- double u = itsUVW[0];
- double v = itsUVW[1];
- itsUVW.backward(1, q);
-
- itsUVW.forward(1, p);
- u -= itsUVW[0];
- v -= itsUVW[1];
- itsUVW.backward(1, p);
-
- // Rotate (u, v) by the position angle and scale with the major
- // and minor axis lengths (FWHM in rad).
- const double uPrime = uScale * (u * cosPhi - v * sinPhi);
- const double vPrime = vScale * (u * sinPhi + v * cosPhi);
-
- // Compute uPrime^2 + vPrime^2 and pre-multiply with -2.0 * PI^2
- // / C^2.
- const double uvPrime = (-2.0 * casa::C::pi * casa::C::pi)
- * (uPrime * uPrime + vPrime * vPrime);
-
- const dcomplex *shiftP = &(itsShiftBuffer[p * itsChannelCount]);
- const dcomplex *shiftQ = &(itsShiftBuffer[q * itsChannelCount]);
- const dcomplex *spectrum = &(itsSpectrumBuffer[0]);
-
- for(size_t ch = 0; ch < itsChannelCount; ++ch)
- {
- // Compute baseline phase shift.
- dcomplex blShift = (*shiftQ) * conj(*shiftP);
- ++shiftP;
- ++shiftQ;
-
- const double ampl = exp(((*itsFreq) * (*itsFreq))
- / (casa::C::c * casa::C::c) * uvPrime);
- ++itsFreq;
-
- blShift *= ampl;
-
- // Compute visibilities.
- *itsBuffer += blShift * (*spectrum);
- ++itsBuffer;
- ++spectrum;
- *itsBuffer += blShift * (*spectrum);
- ++itsBuffer;
- ++spectrum;
- *itsBuffer += blShift * (*spectrum);
- ++itsBuffer;
- ++spectrum;
- *itsBuffer += blShift * (*spectrum);
- ++itsBuffer;
- ++spectrum;
-
- // Move to next channel.
- itsBuffer -= 4;
- itsBuffer.forward(1);
- } // Channels.
- itsFreq -= itsChannelCount;
- itsBuffer.backward(1, itsChannelCount);
- }
-
- // Move to next baseline.
- itsBuffer.forward(2);
- ++itsBaselines;
- } // Baselines.
-
- itsBuffer.backward(2, itsBaselineCount);
- itsBaselines -= itsBaselineCount;
- }
-
- virtual void visit(const Patch&)
- {
- }
-
-private:
- Position itsReference;
- size_t itsStationCount;
- size_t itsBaselineCount;
- size_t itsChannelCount;
- const_cursor<Baseline> itsBaselines;
- const_cursor<double> itsFreq;
- const_cursor<double> itsUVW;
- cursor<dcomplex> itsBuffer;
- vector<dcomplex> itsShiftBuffer;
- vector<dcomplex> itsSpectrumBuffer;
-};
-
void simulate(const Position &reference, const Patch::ConstPtr &patch,
size_t nStation, size_t nBaseline, size_t nChannel,
const_cursor<Baseline> baselines, const_cursor<double> freq,
const_cursor<double> uvw, cursor<dcomplex> vis)
{
- SimulateVisitor visitor(reference, nStation, nBaseline, nChannel, baselines,
+ Simulator simulator(reference, nStation, nBaseline, nChannel, baselines,
freq, uvw, vis);
- visitor.simulate(patch);
-}
-
-namespace
-{
-inline void radec2lmn(const Position &reference, const Position &position,
- cursor<double> lmn)
-{
- const double dRA = position[0] - reference[0];
- const double pDEC = position[1];
- const double rDEC = reference[1];
- const double cDEC = cos(pDEC);
-
- const double l = cDEC * sin(dRA);
- const double m = sin(pDEC) * cos(rDEC) - cDEC * sin(rDEC) * cos(dRA);
-
- lmn[0] = l;
- lmn[1] = m;
- lmn[2] = sqrt(1.0 - l * l - m * m);
-}
-
-inline void phases(size_t nStation, size_t nChannel, const_cursor<double> lmn,
- const_cursor<double> uvw, const_cursor<double> freq, cursor<dcomplex> shift)
-{
- // Compute station phase shifts.
- for(size_t st = 0; st < nStation; ++st)
+ for(size_t i = 0; i < patch->nComponents(); ++i)
{
- const double phase = casa::C::_2pi * (uvw[0] * lmn[0]
- + uvw[1] * lmn[1] + uvw[2] * (lmn[2] - 1.0));
- uvw.forward(1);
-
- for(size_t ch = 0; ch < nChannel; ++ch)
- {
- const double chPhase = phase * (*freq) / casa::C::c;
- *shift = dcomplex(cos(chPhase), sin(chPhase));
- ++freq;
- ++shift;
- } // Channels.
- freq -= nChannel;
- } // Stations.
-}
-
-inline void spectrum(const PointSource &source, size_t nChannel,
- const_cursor<double> freq, cursor<dcomplex> spectrum)
-{
- // Compute source spectrum.
- for(size_t ch = 0; ch < nChannel; ++ch)
- {
- Stokes stokes = source.stokes(*freq);
- ++freq;
-
- spectrum[0] = dcomplex(stokes.I + stokes.Q, 0.0);
- spectrum[1] = dcomplex(stokes.U, stokes.V);
- spectrum[2] = dcomplex(stokes.U, -stokes.V);
- spectrum[3] = dcomplex(stokes.I - stokes.Q, 0.0);
- spectrum += 4;
+ simulator.simulate(patch->component(i));
}
}
-} // Unnamed namespace.
} //# namespace DPPP
} //# namespace LOFAR
diff --git a/CEP/DP3/DPPP/src/Simulator.cc b/CEP/DP3/DPPP/src/Simulator.cc
new file mode 100644
index 00000000000..1e00ed20bca
--- /dev/null
+++ b/CEP/DP3/DPPP/src/Simulator.cc
@@ -0,0 +1,304 @@
+//# Simulator.cc: Compute visibilities for different model components types
+//# (implementation of ModelComponentVisitor).
+//#
+//# Copyright (C) 2012
+//# ASTRON (Netherlands Institute for Radio Astronomy)
+//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//#
+//# This file is part of the LOFAR software suite.
+//# The LOFAR software suite 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 3 of the License, or
+//# (at your option) any later version.
+//#
+//# The LOFAR software suite 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 the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//#
+//# $Id$
+
+#include <lofar_config.h>
+#include <DPPP/Simulator.h>
+#include <DPPP/GaussianSource.h>
+#include <DPPP/PointSource.h>
+#include <casa/BasicSL/Constants.h>
+
+namespace LOFAR
+{
+namespace DPPP
+{
+
+namespace
+{
+// Compute LMN coordinates of \p position relative to \p reference.
+//
+// \param[in] reference
+// Reference position on the celestial sphere.
+// \param[in] position
+// Position of interest on the celestial sphere.
+// \param[in] lmn
+// Pointer to a buffer of (at least) length three into which the computed LMN
+// coordinates will be written.
+void radec2lmn(const Position &reference, const Position &position,
+ cursor<double> lmn);
+
+void phases(size_t nStation, size_t nChannel, const_cursor<double> lmn,
+ const_cursor<double> uvw, const_cursor<double> freq,
+ cursor<dcomplex> shift);
+
+void spectrum(const PointSource &component, size_t nChannel,
+ const_cursor<double> freq, cursor<dcomplex> spectrum);
+} // Unnamed namespace.
+
+Simulator::Simulator(const Position &reference, size_t nStation,
+ size_t nBaseline, size_t nChannel, const_cursor<Baseline> baselines,
+ const_cursor<double> freq, const_cursor<double> uvw,
+ cursor<dcomplex> buffer)
+ : itsReference(reference),
+ itsNStation(nStation),
+ itsNBaseline(nBaseline),
+ itsNChannel(nChannel),
+ itsBaselines(baselines),
+ itsFreq(freq),
+ itsUVW(uvw),
+ itsBuffer(buffer),
+ itsShiftBuffer(nStation * nChannel),
+ itsSpectrumBuffer(nChannel * 4)
+{
+}
+
+void Simulator::simulate(const ModelComponent::ConstPtr &component)
+{
+ component->accept(*this);
+}
+
+void Simulator::visit(const PointSource &component)
+{
+ // Compute LMN coordinates.
+ double lmn[3];
+ radec2lmn(itsReference, component.position(), cursor<double>(lmn));
+
+ // Compute station phase shifts.
+ phases(itsNStation, itsNChannel, const_cursor<double>(lmn),
+ itsUVW, itsFreq, cursor<dcomplex>(&itsShiftBuffer[0]));
+
+ // Compute component spectrum.
+ spectrum(component, itsNChannel, itsFreq,
+ cursor<dcomplex>(&itsSpectrumBuffer[0]));
+
+ // Compute visibilities.
+ for(size_t bl = 0; bl < itsNBaseline; ++bl)
+ {
+ const size_t p = itsBaselines->first;
+ const size_t q = itsBaselines->second;
+
+ if(p != q)
+ {
+ const dcomplex *shiftP = &(itsShiftBuffer[p * itsNChannel]);
+ const dcomplex *shiftQ = &(itsShiftBuffer[q * itsNChannel]);
+ const dcomplex *spectrum = &(itsSpectrumBuffer[0]);
+ for(size_t ch = 0; ch < itsNChannel; ++ch)
+ {
+ // Compute baseline phase shift.
+ const dcomplex blShift = (*shiftQ) * conj(*shiftP);
+ ++shiftP;
+ ++shiftQ;
+
+ // Compute visibilities.
+ *itsBuffer += blShift * (*spectrum);
+ ++itsBuffer;
+ ++spectrum;
+ *itsBuffer += blShift * (*spectrum);
+ ++itsBuffer;
+ ++spectrum;
+ *itsBuffer += blShift * (*spectrum);
+ ++itsBuffer;
+ ++spectrum;
+ *itsBuffer += blShift * (*spectrum);
+ ++itsBuffer;
+ ++spectrum;
+
+ // Move to next channel.
+ itsBuffer -= 4;
+ itsBuffer.forward(1);
+ } // Channels.
+
+ itsBuffer.backward(1, itsNChannel);
+ }
+
+ // Move to next baseline.
+ itsBuffer.forward(2);
+ ++itsBaselines;
+ } // Baselines.
+
+ itsBuffer.backward(2, itsNBaseline);
+ itsBaselines -= itsNBaseline;
+}
+
+void Simulator::visit(const GaussianSource &component)
+{
+ // Compute LMN coordinates.
+ double lmn[3];
+ radec2lmn(itsReference, component.position(), cursor<double>(lmn));
+
+ // Compute station phase shifts.
+ phases(itsNStation, itsNChannel, const_cursor<double>(lmn),
+ itsUVW, itsFreq, cursor<dcomplex>(&itsShiftBuffer[0]));
+
+ // Compute component spectrum.
+ spectrum(component, itsNChannel, itsFreq,
+ cursor<dcomplex>(&itsSpectrumBuffer[0]));
+
+ // Convert position angle from North over East to the angle used to
+ // rotate the right-handed UV-plane.
+ // TODO: Can probably optimize by changing the rotation matrix instead.
+ const double phi = casa::C::pi_2 - component.positionAngle();
+ const double cosPhi = cos(phi);
+ const double sinPhi = sin(phi);
+
+ // Take care of the conversion of axis lengths from FWHM in radians to
+ // sigma.
+ // TODO: Shouldn't the projection from the celestial sphere to the
+ // UV-plane be taken into account here?
+ const double fwhm2sigma = 1.0 / (2.0 * std::sqrt(2.0 * std::log(2.0)));
+ const double uScale = component.majorAxis() * fwhm2sigma;
+ const double vScale = component.minorAxis() * fwhm2sigma;
+
+ for(size_t bl = 0; bl < itsNBaseline; ++bl)
+ {
+ const size_t p = itsBaselines->first;
+ const size_t q = itsBaselines->second;
+
+ if(p != q)
+ {
+ itsUVW.forward(1, q);
+ double u = itsUVW[0];
+ double v = itsUVW[1];
+ itsUVW.backward(1, q);
+
+ itsUVW.forward(1, p);
+ u -= itsUVW[0];
+ v -= itsUVW[1];
+ itsUVW.backward(1, p);
+
+ // Rotate (u, v) by the position angle and scale with the major
+ // and minor axis lengths (FWHM in rad).
+ const double uPrime = uScale * (u * cosPhi - v * sinPhi);
+ const double vPrime = vScale * (u * sinPhi + v * cosPhi);
+
+ // Compute uPrime^2 + vPrime^2 and pre-multiply with -2.0 * PI^2
+ // / C^2.
+ const double uvPrime = (-2.0 * casa::C::pi * casa::C::pi)
+ * (uPrime * uPrime + vPrime * vPrime);
+
+ const dcomplex *shiftP = &(itsShiftBuffer[p * itsNChannel]);
+ const dcomplex *shiftQ = &(itsShiftBuffer[q * itsNChannel]);
+ const dcomplex *spectrum = &(itsSpectrumBuffer[0]);
+
+ for(size_t ch = 0; ch < itsNChannel; ++ch)
+ {
+ // Compute baseline phase shift.
+ dcomplex blShift = (*shiftQ) * conj(*shiftP);
+ ++shiftP;
+ ++shiftQ;
+
+ const double ampl = exp(((*itsFreq) * (*itsFreq))
+ / (casa::C::c * casa::C::c) * uvPrime);
+ ++itsFreq;
+
+ blShift *= ampl;
+
+ // Compute visibilities.
+ *itsBuffer += blShift * (*spectrum);
+ ++itsBuffer;
+ ++spectrum;
+ *itsBuffer += blShift * (*spectrum);
+ ++itsBuffer;
+ ++spectrum;
+ *itsBuffer += blShift * (*spectrum);
+ ++itsBuffer;
+ ++spectrum;
+ *itsBuffer += blShift * (*spectrum);
+ ++itsBuffer;
+ ++spectrum;
+
+ // Move to next channel.
+ itsBuffer -= 4;
+ itsBuffer.forward(1);
+ } // Channels.
+ itsFreq -= itsNChannel;
+ itsBuffer.backward(1, itsNChannel);
+ }
+
+ // Move to next baseline.
+ itsBuffer.forward(2);
+ ++itsBaselines;
+ } // Baselines.
+
+ itsBuffer.backward(2, itsNBaseline);
+ itsBaselines -= itsNBaseline;
+}
+
+namespace
+{
+inline void radec2lmn(const Position &reference, const Position &position,
+ cursor<double> lmn)
+{
+ const double dRA = position[0] - reference[0];
+ const double pDEC = position[1];
+ const double rDEC = reference[1];
+ const double cDEC = cos(pDEC);
+
+ const double l = cDEC * sin(dRA);
+ const double m = sin(pDEC) * cos(rDEC) - cDEC * sin(rDEC) * cos(dRA);
+
+ lmn[0] = l;
+ lmn[1] = m;
+ lmn[2] = sqrt(1.0 - l * l - m * m);
+}
+
+inline void phases(size_t nStation, size_t nChannel, const_cursor<double> lmn,
+ const_cursor<double> uvw, const_cursor<double> freq, cursor<dcomplex> shift)
+{
+ // Compute station phase shifts.
+ for(size_t st = 0; st < nStation; ++st)
+ {
+ const double phase = casa::C::_2pi * (uvw[0] * lmn[0]
+ + uvw[1] * lmn[1] + uvw[2] * (lmn[2] - 1.0));
+ uvw.forward(1);
+
+ for(size_t ch = 0; ch < nChannel; ++ch)
+ {
+ const double chPhase = phase * (*freq) / casa::C::c;
+ *shift = dcomplex(cos(chPhase), sin(chPhase));
+ ++freq;
+ ++shift;
+ } // Channels.
+ freq -= nChannel;
+ } // Stations.
+}
+
+inline void spectrum(const PointSource &component, size_t nChannel,
+ const_cursor<double> freq, cursor<dcomplex> spectrum)
+{
+ // Compute component spectrum.
+ for(size_t ch = 0; ch < nChannel; ++ch)
+ {
+ Stokes stokes = component.stokes(*freq);
+ ++freq;
+
+ spectrum[0] = dcomplex(stokes.I + stokes.Q, 0.0);
+ spectrum[1] = dcomplex(stokes.U, stokes.V);
+ spectrum[2] = dcomplex(stokes.U, -stokes.V);
+ spectrum[3] = dcomplex(stokes.I - stokes.Q, 0.0);
+ spectrum += 4;
+ }
+}
+} // Unnamed namespace.
+
+} //# namespace DPPP
+} //# namespace LOFAR
diff --git a/CEP/DP3/DPPP/src/SourceDBUtil.cc b/CEP/DP3/DPPP/src/SourceDBUtil.cc
index 9217ca7d773..3dcd5c6037a 100644
--- a/CEP/DP3/DPPP/src/SourceDBUtil.cc
+++ b/CEP/DP3/DPPP/src/SourceDBUtil.cc
@@ -110,9 +110,8 @@ vector<Patch::ConstPtr> makePatches(SourceDB &sourceDB,
// Fetch rotation measure attributes (if applicable).
if(src.getInfo().getUseRotationMeasure())
{
- source->setPolarizedFraction(src.getPolarizedFraction());
- source->setPolarizationAngle(src.getPolarizationAngle());
- source->setRotationMeasure(src.getRotationMeasure());
+ source->setRotationMeasure(src.getPolarizedFraction(),
+ src.getPolarizationAngle(), src.getRotationMeasure());
}
componentsList[i].push_back(source);
diff --git a/CEP/DP3/DPPP/src/SubtractMixed.cc b/CEP/DP3/DPPP/src/SubtractMixed.cc
index 54f24daa75f..3359f4d0723 100644
--- a/CEP/DP3/DPPP/src/SubtractMixed.cc
+++ b/CEP/DP3/DPPP/src/SubtractMixed.cc
@@ -42,7 +42,7 @@ void subtract(size_t nBaseline, size_t nChannel,
{
for(size_t ch = 0; ch < nChannel; ++ch)
{
- // Update visibilities.
+ // Subtract weighted model from data.
*data -= static_cast<fcomplex>((*weight) * (*model));
++weight;
++model;
@@ -60,7 +60,7 @@ void subtract(size_t nBaseline, size_t nChannel,
++model;
++data;
- // Move to next channel.
+ // Move to the next channel.
weight -= 4;
weight.forward(1);
model -= 4;
@@ -74,6 +74,7 @@ void subtract(size_t nBaseline, size_t nChannel,
data.backward(1, nChannel);
}
+ // Move to the next baseline.
weight.forward(2);
model.forward(2);
data.forward(2);
--
GitLab