Skip to content
Snippets Groups Projects
Commit 047e805b authored by Joris van Zwieten's avatar Joris van Zwieten
Browse files

Task #3223: Cleaned up Demixer.cc, removed boost::multi_array dependency,... 

Task #3223: Cleaned up Demixer.cc, removed boost::multi_array dependency, several other source files cleaned up and documented.
parent 157f2243
Branches
Tags
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
Showing
with 1157 additions and 1455 deletions
.gitattributes
+ 2
− 0
View file @ 047e805b
......@@ -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
......
CEP/DP3/DPPP/include/DPPP/BBSExpr.h deleted 100644 → 0
+ 0
− 111
View file @ 157f2243
//# 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
CEP/DP3/DPPP/include/DPPP/CMakeLists.txt
+ 1
− 1
View file @ 047e805b
......@@ -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)
......
CEP/DP3/DPPP/include/DPPP/Demixer.h
+ 20
− 39
View file @ 047e805b
......@@ -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
......
CEP/DP3/DPPP/include/DPPP/EstimateMixed.h
+ 1
− 1
View file @ 047e805b
......@@ -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
......
CEP/DP3/DPPP/include/DPPP/ModelComponentVisitor.h
+ 0
− 2
View file @ 047e805b
......@@ -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;
};
// @}
......
CEP/DP3/DPPP/include/DPPP/Patch.h
+ 15
− 3
View file @ 047e805b
......@@ -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
......
CEP/DP3/DPPP/include/DPPP/PointSource.h
+ 2
− 7
View file @ 047e805b
......@@ -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;
};
// @}
......
CEP/DP3/DPPP/include/DPPP/Simulator.h 0 → 100644
+ 75
− 0
View file @ 047e805b
//# 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
CEP/DP3/DPPP/src/Apply.cc
+ 37
− 48
View file @ 047e805b
......@@ -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.
......
CEP/DP3/DPPP/src/BBSExpr.cc deleted 100644 → 0
+ 0
− 221
View file @ 157f2243
//# 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
CEP/DP3/DPPP/src/CMakeLists.txt
+ 1
− 1
View file @ 047e805b
......@@ -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
......
CEP/DP3/DPPP/src/Demixer.cc
+ 357
− 447
View file @ 047e805b
This diff is collapsed.
CEP/DP3/DPPP/src/EstimateMixed.cc
+ 145
− 179
View file @ 047e805b
......@@ -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;
......
CEP/DP3/DPPP/src/Patch.cc
+ 16
− 27
View file @ 047e805b
......@@ -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,11 +44,8 @@ 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)
{
......@@ -75,6 +63,7 @@ void Patch::computePosition()
itsPosition[0] = atan2(y, x);
itsPosition[1] = asin(z);
}
}
} //# namespace DPPP
} //# namespace LOFAR
CEP/DP3/DPPP/src/PhaseShift.cc
+ 24
− 24
View file @ 047e805b
CEP/DP3/DPPP/src/PointSource.cc
+ 8
− 13
View file @ 047e805b
......@@ -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
......
CEP/DP3/DPPP/src/Simulate.cc
+ 147
− 328
View file @ 047e805b
......@@ -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,40 +35,112 @@ namespace LOFAR
namespace DPPP
{
namespace
void splitUVW(size_t nStation, size_t nBaseline,
const_cursor<Baseline> baselines, const_cursor<double> uvw,
cursor<double> split)
{
// 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);
// 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);
void phases(size_t nStation, size_t nChannel, const_cursor<double> lmn,
const_cursor<double> uvw, const_cursor<double> freq,
cursor<dcomplex> shift);
// Flags to keep track of the stations for which the UVW coordinates are
// known.
vector<bool> flag(nStation, true);
void spectrum(const PointSource &source, size_t nChannel,
const_cursor<double> freq, cursor<dcomplex> spectrum);
} // Unnamed namespace.
// 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 = tmp->first;
const size_t q = tmp->second;
void splitUVW(size_t nStation, size_t nBaseline,
const_cursor<Baseline> baselines, const_cursor<double> uvw,
cursor<double> split)
if(p != q)
{
if(flag[p])
{
flag[p] = false;
++nReachable;
}
if(flag[q])
{
flag[q] = false;
++nReachable;
}
}
if(nReachable == nStation)
{
break;
}
// Move to next baseline.
++tmp;
}
// 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);
vector<bool> flag(nStation, false);
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;
flag[0] = true;
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;
......@@ -101,6 +169,8 @@ void splitUVW(size_t nStation, size_t nBaseline,
known.backward(1, q);
flag[p] = true;
}
++nFound;
}
// Move to next baseline.
......@@ -108,9 +178,22 @@ void splitUVW(size_t nStation, size_t nBaseline,
++baselines;
} // Baselines.
ASSERTSTR(static_cast<size_t>(std::count(flag.begin(), flag.end(), true))
== flag.size(), "Unable to split baseline UVW coordinates into station"
" UVW coordinates.");
// 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)
for(size_t i = 0; i < patch->nComponents(); ++i)
{
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 &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
CEP/DP3/DPPP/src/Simulator.cc 0 → 100644
+ 304
− 0
View file @ 047e805b
//# 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
CEP/DP3/DPPP/src/SourceDBUtil.cc
+ 2
− 3
View file @ 047e805b
......@@ -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);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment