diff --git a/CEP/DP3/DPPP/src/SolTab.cc b/CEP/DP3/DPPP/src/SolTab.cc
index 54468913db0a69fedfdf9dc581066bc56603d9a7..625e2d7c580dbd0ca0fafaa23d97a44dfe308256 100644
--- a/CEP/DP3/DPPP/src/SolTab.cc
+++ b/CEP/DP3/DPPP/src/SolTab.cc
@@ -125,7 +125,12 @@ namespace LOFAR {
}
// Add weights
- H5::DataSet weightset = createDataSet("weight", H5::PredType::IEEE_F64LE,
+ hid_t halffloat = H5Tcopy(H5T_IEEE_F32BE);
+ H5Tset_fields(halffloat, 15, 10, 5, 0, 10);
+ H5Tset_size(halffloat, 2);
+ H5Tset_ebias(halffloat, 15);
+ H5Tlock(halffloat);
+ H5::DataSet weightset = createDataSet("weight", H5::DataType(halffloat),
dataspace);
// If weights are empty, write ones everywhere
diff --git a/CEP/DP3/DPPP/test/tH5Parm.cc b/CEP/DP3/DPPP/test/tH5Parm.cc
index 4291800d1dbc595fa0d7ed8a33ef3d6c2d103440..ddec947dc907e57c759c0fb2b7eb1d7ce1586bdb 100644
--- a/CEP/DP3/DPPP/test/tH5Parm.cc
+++ b/CEP/DP3/DPPP/test/tH5Parm.cc
@@ -9,7 +9,7 @@
using namespace std;
using namespace LOFAR;
-void checkAxes(H5Parm::SolTab& soltab) {
+void checkAxes(H5Parm::SolTab& soltab, size_t ntimes) {
ASSERT(soltab.nAxes()==3);
ASSERT(soltab.hasAxis("ant"));
ASSERT(soltab.hasAxis("time"));
@@ -18,12 +18,13 @@ void checkAxes(H5Parm::SolTab& soltab) {
ASSERT(soltab.getAxis(1).name=="time");
ASSERT(soltab.getAxis(2).name=="bla");
ASSERT(soltab.getAxis(0).size==3);
- ASSERT(soltab.getAxis(1).size==4);
+ ASSERT(soltab.getAxis(1).size==ntimes);
ASSERT(soltab.getAxis(2).size==1);
}
int main(int, char**) {
{
+ size_t ntimes=42;
{
// Create a new H5Parm
cout<<"Create tH5Parm_tmp.h5"<<endl;
@@ -49,7 +50,7 @@ int main(int, char**) {
vector<H5Parm::AxisInfo> axes;
axes.push_back(H5Parm::AxisInfo("ant",3));
- axes.push_back(H5Parm::AxisInfo("time",4));
+ axes.push_back(H5Parm::AxisInfo("time",ntimes));
axes.push_back(H5Parm::AxisInfo("bla",1));
cout<<"Create new SolTab"<<endl;
@@ -62,14 +63,19 @@ int main(int, char**) {
// Check the axes
H5Parm::SolTab soltab = h5parm.getSolTab("mysol");
ASSERT(soltab.getType()=="mytype");
- checkAxes(soltab);
+ checkAxes(soltab, ntimes);
// Add some data
- vector<double> vals(3*4);
- for (size_t ant=0; ant<3; ++ant)
- for (size_t time=0; time<4; ++time)
- vals[ant*4+time]=10*ant+time;
- soltab.setValues(vals, vector<double>(), "CREATE with DPPP");
+ vector<double> vals(3*ntimes);
+ vector<double> weights(3*ntimes);
+ for (size_t ant=0; ant<3; ++ant) {
+ for (size_t time=0; time<ntimes; ++time) {
+ vals[ant*ntimes+time]=10*ant+time;
+ weights[ant*ntimes+time]=0.4;
+ }
+ }
+
+ soltab.setValues(vals, weights, "CREATE with DPPP");
// Add metadata for stations
vector<string> someAntNames;
@@ -80,9 +86,9 @@ int main(int, char**) {
// Add metadata for times
vector<double> times;
- times.push_back(57878.5);
- times.push_back(57880.5);
- times.push_back(57882.5);
+ for (size_t time=0; time<ntimes; ++time) {
+ times.push_back(57878.5+2.0*time);
+ }
soltab.setTimes(times);
// Add metadata for freqs;
@@ -117,13 +123,13 @@ int main(int, char**) {
// Check the axes
H5Parm::SolTab soltab = h5parm.getSolTab("mysol");
ASSERT(soltab.getType()=="mytype");
- checkAxes(soltab);
+ checkAxes(soltab, ntimes);
cout<<"read some data"<<endl;
double starttime = 57878.49999;
hsize_t starttimeindex = soltab.getTimeIndex(starttime);
cout<<"starttimeindex="<<starttimeindex<<endl;
- vector<double> val = soltab.getValues("Antenna2", starttimeindex, 4);
+ vector<double> val = soltab.getValues("Antenna2", starttimeindex, ntimes);
ASSERT(casa::near(val[0],10.));
ASSERT(casa::near(val[1],11.));
ASSERT(casa::near(val[2],12.));
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Constraint.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Constraint.h
index 7f215b28d98a0870fc15b470f1fff31e3dd2caa1..49ff5b891f408f12e6f58ca3bc764233c03d890e 100644
--- a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Constraint.h
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Constraint.h
@@ -21,6 +21,7 @@ public:
{
public:
std::vector<double> vals;
+ std::vector<double> weights;
std::string axes; // Comma-separated string with axis names, fastest varying last
std::vector<size_t> dims;
std::string name;
@@ -65,7 +66,30 @@ public:
std::vector<std::vector<dcomplex> >& solutions,
double time) = 0;
+ /**
+ * Initialize the dimensions for the constraint. Should be overridden when
+ * something more than assigning dimensions is needed (e.g. resizing vectors).
+ * Weights are initialized to 1. here.
+ */
+ virtual void InitializeDimensions(size_t nAntennas,
+ size_t nDirections,
+ size_t nChannelBlocks)
+ {
+ _nAntennas = nAntennas;
+ _nDirections = nDirections;
+ _nChannelBlocks = nChannelBlocks;
+ }
+
+ /**
+ * Set weights. The vector should contain an array of size nAntennas * nChannelBlocks,
+ * where the channel index varies fastest.
+ */
+ virtual void SetWeights(std::vector<double> &) {}
+
virtual void showTimings (std::ostream&, double) const {}
+
+protected:
+ size_t _nAntennas, _nDirections, _nChannelBlocks;
};
/**
@@ -123,11 +147,8 @@ class CoreConstraint : public Constraint
public:
CoreConstraint() { }
- void initialize(size_t nAntennas, size_t nDirections, size_t nChannelBlocks, const std::set<size_t>& coreAntennas)
+ void initialize(const std::set<size_t>& coreAntennas)
{
- _nAntennas = nAntennas;
- _nDirections = nDirections;
- _nChannelBlocks = nChannelBlocks;
_coreAntennas = coreAntennas;
}
@@ -136,7 +157,6 @@ public:
double time);
private:
- size_t _nAntennas, _nDirections, _nChannelBlocks;
std::set<size_t> _coreAntennas;
};
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h
index f4ec891d3eefb3ca5025d77aa4d599a3d6cdd16a..d74e77465bbed6567cc452721f3977e9e14a79dc 100644
--- a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h
@@ -141,10 +141,13 @@ namespace LOFAR {
uint itsSolInt;
uint itsStepInSolInt;
uint itsNChan;
+ vector<size_t> itsChanBlockStart; // For each channel block, the index in the channels at which this channel block starts
vector<double> itsChanBlockFreqs;
vector<vector<string> > itsDirections; // For each direction, a vector of patches
vector<casacore::CountedPtr<Constraint> > itsConstraints;
+ vector<double> itsWeights;
+
UVWFlagger itsUVWFlagStep;
ResultStep::ShPtr itsDataResultStep; // Result step for data after UV-flagging
vector<Predict> itsPredictSteps;
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/RotationAndDiagonalConstraint.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/RotationAndDiagonalConstraint.h
index 996f95130c2860eded0717cdeb908fcceff87f1d..a7caac05a9d7b44b7a22c2448fd927a099bfe9fd 100644
--- a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/RotationAndDiagonalConstraint.h
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/RotationAndDiagonalConstraint.h
@@ -14,16 +14,19 @@ namespace LOFAR {
class RotationAndDiagonalConstraint : public Constraint
{
public:
- RotationAndDiagonalConstraint();
+ RotationAndDiagonalConstraint() {};
virtual std::vector<Result> Apply(
std::vector<std::vector<dcomplex> >& solutions,
double time);
- void initialize(size_t nAntennas, size_t nDirections, size_t nChannelBlocks);
+ virtual void InitializeDimensions(size_t nAntennas,
+ size_t nDirections,
+ size_t nChannelBlocks);
+
+ virtual void SetWeights(std::vector<double>&);
private:
- size_t _nAntennas, _nDirections, _nChannelBlocks;
std::vector<Constraint::Result> _resTemplate;
};
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/RotationConstraint.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/RotationConstraint.h
index ddc25f71a59268d4c19c19654d5d3879d8cefaf3..034b9c7ba4a0aab3ae35c0455f34a0ee1377e959 100644
--- a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/RotationConstraint.h
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/RotationConstraint.h
@@ -14,19 +14,22 @@ namespace LOFAR {
class RotationConstraint : public Constraint
{
public:
- RotationConstraint();
+ RotationConstraint() {};
virtual std::vector<Result> Apply(
std::vector<std::vector<dcomplex> >& solutions,
double time);
- void initialize(size_t nAntennas, size_t nDirections, size_t nChannelBlocks);
+ virtual void InitializeDimensions(size_t nAntennas,
+ size_t nDirections,
+ size_t nChannelBlocks);
+
+ virtual void SetWeights(std::vector<double>&);
/* Compute the rotation from a 2x2 full jones solution */
static double get_rotation(std::complex<double>* data);
private:
- size_t _nAntennas, _nDirections, _nChannelBlocks;
std::vector<Constraint::Result> _resTemplate;
};
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/ScreenConstraint.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/ScreenConstraint.h
index d22abd5353344c902953764d754304f68586207d..1ac33109aa31a3a4eb0089bb4c9c2271453d5b37 100644
--- a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/ScreenConstraint.h
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/ScreenConstraint.h
@@ -21,7 +21,11 @@ class ScreenConstraint : public Constraint
public:
ScreenConstraint(const ParameterSet& parset,
const string& prefix);
- void initialize(size_t nAntennas, size_t nDirections, size_t nChannelBlocks, const double* frequencies);
+
+ /** Initialize metadata with frequencies, resize some members.
+ * Should be called after InitializeDimensions.
+ */
+ void initialize(const double* frequencies);
virtual std::vector<Constraint::Result> Apply(std::vector<std::vector<MultiDirSolver::DComplex> >& solutions,double time);
virtual void CalculatePiercepoints();
@@ -41,7 +45,6 @@ public:
}
void getPPValue(std::vector<std::vector<std::complex<double> > >&, size_t, size_t, double&,double&) const;
private:
- size_t _nAntennas, _nDirections, _nChannelBlocks;
std::vector<std::vector<double> > itsAntennaPos;
std::vector<std::vector<double> > itsSourcePos;
std::vector<double> itsFrequencies;
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/TECConstraint.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/TECConstraint.h
index 647935622d96c3436e68d0bd1a14f6fdb469eb14..342ede6312fb738dcb2a456e59e063276a863cf3 100644
--- a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/TECConstraint.h
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/TECConstraint.h
@@ -1,6 +1,8 @@
#ifndef TEC_CONSTRAINT_H
#define TEC_CONSTRAINT_H
+#include <vector>
+
#ifdef AOPROJECT
#include "Constraint.h"
#include "PhaseFitter.h"
@@ -25,16 +27,20 @@ public:
TECConstraintBase(Mode mode);
- void initialize(size_t nAntennas, size_t nDirections,
- size_t nChannelBlocks, const double* frequencies);
+ /** Initialize metadata with frequencies, resize some members.
+ * Should be called after InitializeDimensions.
+ */
+ void initialize(const double* frequencies);
+ /** Propagate weights to the phase fitters */
+ virtual void SetWeights(std::vector<double>& weights);
+
protected:
virtual void initializeChild() { }
void applyReferenceAntenna(std::vector<std::vector<dcomplex> >& solutions) const;
Mode _mode;
- size_t _nAntennas, _nDirections, _nChannelBlocks;
std::vector<PhaseFitter> _phaseFitters;
};
diff --git a/CEP/DP3/DPPP_DDECal/src/DDECal.cc b/CEP/DP3/DPPP_DDECal/src/DDECal.cc
index 6e19eb7325a2a0f922e10f868ee52420813bc342..c1329eb0bc13522639fc524c203270b031a19b0b 100644
--- a/CEP/DP3/DPPP_DDECal/src/DDECal.cc
+++ b/CEP/DP3/DPPP_DDECal/src/DDECal.cc
@@ -327,6 +327,7 @@ namespace LOFAR {
itsConstraintSols.resize(nSolTimes);
size_t nChannelBlocks = info().nchan()/itsNChan;
+ itsChanBlockStart.resize(nChannelBlocks+1);
itsChanBlockFreqs.resize(nChannelBlocks);
for(size_t chBlock=0; chBlock!=nChannelBlocks; ++chBlock) {
const size_t
@@ -337,19 +338,29 @@ namespace LOFAR {
info().chanFreqs().data()+channelIndexStart,
info().chanFreqs().data()+channelIndexEnd,
0.0) / curChannelBlockSize;
+ itsChanBlockStart[chBlock] = channelIndexStart;
itsChanBlockFreqs[chBlock] = meanfreq;
}
+ itsChanBlockStart[itsChanBlockStart.size()-1] = info().nchan();
+
+ itsWeights.resize(itsChanBlockFreqs.size()*info().nantenna());
for (uint i=0; i<itsConstraints.size();++i) {
+ // Initialize the constraint with some common metadata
+ itsConstraints[i]->InitializeDimensions(info().antennaNames().size(),
+ itsDirections.size(),
+ nChannelBlocks);
+
// Different constraints need different information. Determine if the constraint is
// of a type that needs more information, and if so initialize the constraint.
CoreConstraint* coreConstraint = dynamic_cast<CoreConstraint*>(itsConstraints[i].get());
if(coreConstraint != 0)
{
+ // Take antenna with index 0 as reference station
double
- refX = antennaPos[i][0],
- refY = antennaPos[i][1],
- refZ = antennaPos[i][2];
+ refX = antennaPos[0][0],
+ refY = antennaPos[0][1],
+ refZ = antennaPos[0][2];
std::set<size_t> coreAntennaIndices;
const double coreDistSq = itsCoreConstraint*itsCoreConstraint;
for(size_t ant=0; ant!=antennaPos.size(); ++ant)
@@ -362,21 +373,13 @@ namespace LOFAR {
if(distSq <= coreDistSq)
coreAntennaIndices.insert(ant);
}
- coreConstraint->initialize(info().antennaNames().size(),
- itsDirections.size(),
- info().nchan(),
- coreAntennaIndices);
+ coreConstraint->initialize(coreAntennaIndices);
}
ScreenConstraint* screenConstraint = dynamic_cast<ScreenConstraint*>(itsConstraints[i].get());
if(screenConstraint != 0)
{
- screenConstraint->initialize(
- info().antennaNames().size(),
- itsDirections.size(),
- nChannelBlocks,
- &(itsChanBlockFreqs[0])
- );
+ screenConstraint->initialize(&(itsChanBlockFreqs[0]));
screenConstraint->setAntennaPositions(antennaPos);
screenConstraint->setDirections(sourcePositions);
screenConstraint->initPiercePoints();
@@ -406,21 +409,7 @@ namespace LOFAR {
TECConstraintBase* tecConstraint = dynamic_cast<TECConstraintBase*>(itsConstraints[i].get());
if(tecConstraint != 0)
{
- tecConstraint->initialize(info().antennaNames().size(),
- itsDirections.size(),
- nChannelBlocks,
- &(itsChanBlockFreqs[0]));
- }
-
- RotationAndDiagonalConstraint* rotationAndDiagonalConstraint = dynamic_cast<RotationAndDiagonalConstraint*>(itsConstraints[i].get());
- if(rotationAndDiagonalConstraint != 0)
- {
- rotationAndDiagonalConstraint->initialize(info().antennaNames().size(), itsDirections.size(), nChannelBlocks);
- }
- RotationConstraint* rotationConstraint = dynamic_cast<RotationConstraint*>(itsConstraints[i].get());
- if(rotationConstraint != 0)
- {
- rotationConstraint->initialize(info().antennaNames().size(), itsDirections.size(), nChannelBlocks);
+ tecConstraint->initialize(&(itsChanBlockFreqs[0]));
}
}
@@ -600,7 +589,15 @@ namespace LOFAR {
const size_t nCh = info().nchan();
const size_t nCr = 4;
+ size_t nchanblocks = itsChanBlockFreqs.size();
+ size_t chanblock = 0;
+
+ double weightFactor = 1./(nCh*(info().nantenna()-1)*nCr*itsSolInt);
+
for (size_t ch=0; ch<nCh; ++ch) {
+ if (ch == itsChanBlockStart[chanblock+1]) {
+ chanblock++;
+ }
for (size_t bl=0; bl<nBl; ++bl) {
for (size_t cr=0; cr<nCr; ++cr) {
if (itsBufs[itsStepInSolInt].getFlags().data()[bl*nCr*nCh+ch*nCr+cr]) {
@@ -611,26 +608,39 @@ namespace LOFAR {
}
} else {
// Premultiply non-flagged data with sqrt(weight)
- double weight = sqrt(itsBufs[itsStepInSolInt].getWeights().data()[bl*nCr*nCh+ch*nCr+cr]);
- itsDataPtrs[itsStepInSolInt][bl*nCr*nCh+ch*nCr+cr] *= weight;
+ double weight = itsBufs[itsStepInSolInt].getWeights().data()[bl*nCr*nCh+ch*nCr+cr];
+ itsDataPtrs[itsStepInSolInt][bl*nCr*nCh+ch*nCr+cr] *= sqrt(weight);
+ itsWeights[info().getAnt1()[bl]*nchanblocks + chanblock] += weight;
+ itsWeights[info().getAnt2()[bl]*nchanblocks + chanblock] += weight;
for (size_t dir=0; dir<itsPredictSteps.size(); ++dir) {
- itsModelDataPtrs[itsStepInSolInt][dir][bl*nCr*nCh+ch*nCr+cr] *= weight;
+ itsModelDataPtrs[itsStepInSolInt][dir][bl*nCr*nCh+ch*nCr+cr] *= sqrt(weight);
}
}
}
}
}
+ for (auto& weight: itsWeights) {
+ weight *= weightFactor;
+ }
+
itsTimerPredict.stop();
itsAvgTime += itsAvgTime + bufin.getTime();
if (itsStepInSolInt==itsSolInt-1) {
+ for (uint constraint_num = 0; constraint_num < itsConstraints.size(); ++constraint_num) {
+ itsConstraints[constraint_num]->SetWeights(itsWeights);
+ }
+
doSolve();
- itsStepInSolInt=0;
+
+ // Clean up, prepare for next iteration
+ itsStepInSolInt=0;
itsAvgTime=0;
for (size_t dir=0; dir<itsResultSteps.size(); ++dir) {
itsResultSteps[dir]->clear();
+ itsWeights.assign(itsWeights.size(), 0.);
}
} else {
itsStepInSolInt++;
@@ -822,9 +832,22 @@ namespace LOFAR {
nextpos);
}
+ // Put solution weights in a contiguous piece of memory
+ vector<double> weights;
+ if (!itsConstraintSols[0][constraintNum][solNameNum].weights.empty()) {
+ weights.resize(numSols);
+ vector<double>::iterator nextpos = weights.begin();
+ for (uint time=0; time<itsSols.size(); ++time) {
+ nextpos = std::copy(
+ itsConstraintSols[time][constraintNum][solNameNum].weights.begin(),
+ itsConstraintSols[time][constraintNum][solNameNum].weights.end(),
+ nextpos);
+ }
+ }
+
string solTabName = firstResult.name+"000";
H5Parm::SolTab soltab = itsH5Parm.createSolTab(solTabName, firstResult.name, axes);
- soltab.setValues(sols, vector<double>(),
+ soltab.setValues(sols, weights,
"CREATE by DPPP\n" +
Version::getInfo<DPPPVersion>("DPPP", "top") + "\n" +
"step " + itsName + " in parset: \n" +
diff --git a/CEP/DP3/DPPP_DDECal/src/RotationAndDiagonalConstraint.cc b/CEP/DP3/DPPP_DDECal/src/RotationAndDiagonalConstraint.cc
index b42267e6ae7725d506f116fcaf942a27efdd1095..ef9966bf6f2283bcafdd312c5212f174f0a5db90 100644
--- a/CEP/DP3/DPPP_DDECal/src/RotationAndDiagonalConstraint.cc
+++ b/CEP/DP3/DPPP_DDECal/src/RotationAndDiagonalConstraint.cc
@@ -10,19 +10,16 @@ using namespace std;
namespace LOFAR {
-RotationAndDiagonalConstraint::RotationAndDiagonalConstraint():
- _nAntennas(0), _nDirections(0)
-{
-}
+void RotationAndDiagonalConstraint::InitializeDimensions(size_t nAntennas,
+ size_t nDirections,
+ size_t nChannelBlocks) {
+ Constraint::InitializeDimensions(nAntennas, nDirections, nChannelBlocks);
-void RotationAndDiagonalConstraint::initialize(size_t nAntennas, size_t nDirections, size_t nChannelBlocks) {
- _nAntennas = nAntennas;
- _nDirections = nDirections;
- assert(nDirections == 1);
- _nChannelBlocks = nChannelBlocks;
+ assert(_nDirections == 1);
_resTemplate.resize(3);
_resTemplate[0].vals.resize(_nAntennas*_nChannelBlocks);
+ _resTemplate[0].weights.resize(_nAntennas*_nChannelBlocks);
_resTemplate[0].axes="ant,freq";
_resTemplate[0].dims.resize(2);
_resTemplate[0].dims[0]=_nAntennas;
@@ -30,6 +27,7 @@ void RotationAndDiagonalConstraint::initialize(size_t nAntennas, size_t nDirecti
_resTemplate[0].name="rotation";
_resTemplate[1].vals.resize(_nAntennas*_nChannelBlocks*2);
+ _resTemplate[1].weights.resize(_nAntennas*_nChannelBlocks*2);
_resTemplate[1].axes="ant,freq,pol";
_resTemplate[1].dims.resize(3);
_resTemplate[1].dims[0]=_nAntennas;
@@ -41,6 +39,20 @@ void RotationAndDiagonalConstraint::initialize(size_t nAntennas, size_t nDirecti
_resTemplate[2].name="phase";
}
+void RotationAndDiagonalConstraint::SetWeights(vector<double>& weights) {
+ _resTemplate[0].weights.assign(weights.begin(), weights.end());
+
+ // Duplicate weights for two polarizations
+ _resTemplate[1].weights.resize(_nAntennas*_nChannelBlocks*2);
+ size_t indexInWeights = 0;
+ for (auto weight: weights) {
+ _resTemplate[1].weights[indexInWeights++] = weight;
+ _resTemplate[1].weights[indexInWeights++] = weight;
+ }
+
+ _resTemplate[2].weights = _resTemplate[1].weights;
+}
+
vector<Constraint::Result> RotationAndDiagonalConstraint::Apply(
vector<vector<dcomplex> >& solutions, double) {
// Convert to circular
@@ -55,7 +67,7 @@ vector<Constraint::Result> RotationAndDiagonalConstraint::Apply(
double angle = RotationConstraint::get_rotation(data);
// Restrict angle between -pi/2 and pi/2
- // Add 2pi to make sure that fmod doesn't see negative numbers
+ // Add 2pi to make sure that fmod doesn't see negative numbers
angle = fmod(angle + 3.5*M_PI, M_PI) - 0.5*M_PI;
_resTemplate[0].vals[ant*_nChannelBlocks + ch] = angle;
diff --git a/CEP/DP3/DPPP_DDECal/src/RotationConstraint.cc b/CEP/DP3/DPPP_DDECal/src/RotationConstraint.cc
index 2b14c6578e4d05fe80efa9e8a9fb4b24ae51ef35..913a206182591d909417ef0ea840e516c66de48a 100644
--- a/CEP/DP3/DPPP_DDECal/src/RotationConstraint.cc
+++ b/CEP/DP3/DPPP_DDECal/src/RotationConstraint.cc
@@ -9,16 +9,12 @@ using namespace std;
namespace LOFAR {
-RotationConstraint::RotationConstraint():
- _nAntennas(0), _nDirections(0)
-{
-}
+ void RotationConstraint::InitializeDimensions(size_t nAntennas,
+ size_t nDirections,
+ size_t nChannelBlocks) {
+ Constraint::InitializeDimensions(nAntennas, nDirections, nChannelBlocks);
-void RotationConstraint::initialize(size_t nAntennas, size_t nDirections, size_t nChannelBlocks) {
- _nAntennas = nAntennas;
- _nDirections = nDirections;
- assert(nDirections == 1);
- _nChannelBlocks = nChannelBlocks;
+ assert(_nDirections == 1);
_resTemplate.resize(1);
_resTemplate[0].vals.resize(_nAntennas*_nChannelBlocks);
@@ -29,6 +25,10 @@ void RotationConstraint::initialize(size_t nAntennas, size_t nDirections, size_t
_resTemplate[0].name="rotation";
}
+void RotationConstraint::SetWeights(vector<double>& weights) {
+ _resTemplate[0].weights.assign(weights.begin(), weights.end());
+}
+
double RotationConstraint::get_rotation(std::complex<double>* data) {
// Convert to circular
complex<double> ll, rr;
diff --git a/CEP/DP3/DPPP_DDECal/src/ScreenConstraint.cc b/CEP/DP3/DPPP_DDECal/src/ScreenConstraint.cc
index 5e60d3bcbbce289ddd376ed598549996b972f579..24a78330a65344d0f072c44d5ff2609f88e20213 100644
--- a/CEP/DP3/DPPP_DDECal/src/ScreenConstraint.cc
+++ b/CEP/DP3/DPPP_DDECal/src/ScreenConstraint.cc
@@ -10,10 +10,8 @@ const size_t ScreenConstraint::maxIter=30;
ScreenConstraint::ScreenConstraint(const ParameterSet& parset,
const string& prefix)
:
- _nAntennas(0),
- _nDirections(0),
- _nChannelBlocks(0),
- itsCurrentTime(0)
+ itsCurrentTime(0),
+ itsIter(0)
{
cout<<"=========="<<(prefix + "order")<<"========\n";
itsBeta=parset.getDouble (prefix + "beta", 5./3.);
@@ -25,10 +23,7 @@ ScreenConstraint::ScreenConstraint(const ParameterSet& parset,
itsDebugMode=parset.getInt(prefix + "debug", 0);
}
-void ScreenConstraint::initialize(size_t nAntennas, size_t nDirections, size_t nChannelBlocks, const double* frequencies) {
- _nAntennas = nAntennas;
- _nDirections = nDirections;
- _nChannelBlocks = nChannelBlocks;
+void ScreenConstraint::initialize(const double* frequencies) {
itsFrequencies.resize(_nChannelBlocks);
itsprevsol.assign(_nDirections*_nAntennas,-999.);
std::memcpy( itsFrequencies.data(),frequencies, sizeof(double) * _nChannelBlocks);
diff --git a/CEP/DP3/DPPP_DDECal/src/TECConstraint.cc b/CEP/DP3/DPPP_DDECal/src/TECConstraint.cc
index a7b588667f255be54a2411faa15dff5a32597dc3..fbaf7d5936c10192b06b811f1ffa05240e39e87a 100644
--- a/CEP/DP3/DPPP_DDECal/src/TECConstraint.cc
+++ b/CEP/DP3/DPPP_DDECal/src/TECConstraint.cc
@@ -6,19 +6,15 @@
#include <Common/OpenMP.h>
#endif
+#include <assert.h>
+
TECConstraintBase::TECConstraintBase(Mode mode) :
_mode(mode),
- _nAntennas(0),
- _nDirections(0),
- _nChannelBlocks(0),
_phaseFitters()
{
}
-void TECConstraintBase::initialize(size_t nAntennas, size_t nDirections, size_t nChannelBlocks, const double* frequencies) {
- _nAntennas = nAntennas;
- _nDirections = nDirections;
- _nChannelBlocks = nChannelBlocks;
+void TECConstraintBase::initialize(const double* frequencies) {
_phaseFitters.resize(
#ifdef AOPROJECT
omp_get_max_threads()
@@ -30,13 +26,31 @@ void TECConstraintBase::initialize(size_t nAntennas, size_t nDirections, size_t
for(size_t i=0; i!=_phaseFitters.size(); ++i)
{
_phaseFitters[i].SetChannelCount(_nChannelBlocks);
- std::memcpy(_phaseFitters[i].FrequencyData(), frequencies, sizeof(double) * _nChannelBlocks);
-
- // TODO this should set the weights of the phase fitter!
+ std::memcpy(_phaseFitters[i].FrequencyData(), frequencies,
+ sizeof(double) * _nChannelBlocks);
}
initializeChild();
}
+void TECConstraintBase::SetWeights(std::vector<double> &weights) {
+ std::vector<double> chanBlockWeights(_nChannelBlocks, 0.);
+#ifndef NDEBUG
+ assert(weights.size() == _nAntennas*_nChannelBlocks);
+#endif
+
+ size_t weightspos = 0;
+ for (size_t ant=0; ant < _nAntennas; ++ant) {
+ for (size_t chanBlock=0; chanBlock < _nChannelBlocks; ++chanBlock) {
+ chanBlockWeights[chanBlock] += weights[weightspos++];
+ }
+ }
+
+ for(size_t thread=0; thread!=_phaseFitters.size(); ++thread) {
+ std::copy(chanBlockWeights.begin(), chanBlockWeights.end(),
+ _phaseFitters[thread].WeightData());
+ }
+}
+
void ApproximateTECConstraint::initializeChild()
{
_pwFitters.resize(
diff --git a/CEP/DP3/DPPP_DDECal/test/tRotationConstraint.cc b/CEP/DP3/DPPP_DDECal/test/tRotationConstraint.cc
index a9fb2e35d20563f2c8764416743af8c3136386de..b24ce1b26d8afda49176d4e50bafc08b2f0e5b8f 100644
--- a/CEP/DP3/DPPP_DDECal/test/tRotationConstraint.cc
+++ b/CEP/DP3/DPPP_DDECal/test/tRotationConstraint.cc
@@ -15,7 +15,7 @@ using namespace LOFAR;
void test_rotation() {
RotationConstraint constraint;
- constraint.initialize(1, 1, 1);
+ constraint.InitializeDimensions(1, 1, 1);
vector<vector<complex<double> > > onesolution(1);
onesolution[0].resize(4);
@@ -46,7 +46,7 @@ void test_rotation() {
void test_rotation_and_diagonal() {
RotationAndDiagonalConstraint constraint;
- constraint.initialize(1, 1, 1);
+ constraint.InitializeDimensions(1, 1, 1);
vector<vector<complex<double> > > onesolution(1);
onesolution[0].resize(4);