diff --git a/ddecal/gain_solvers/DiagonalLowRankSolver.cc b/ddecal/gain_solvers/DiagonalLowRankSolver.cc
index 83ef225312824277f429ad74dd5b789ab0c76b7c..c6c499157fdf168a38e8a5d38664dc3e79addc1c 100644
--- a/ddecal/gain_solvers/DiagonalLowRankSolver.cc
+++ b/ddecal/gain_solvers/DiagonalLowRankSolver.cc
@@ -51,7 +51,7 @@ float DominantEigenPairNear(const xt::xtensor<std::complex<float>, 2>& matrix,
}
DiagonalLowRankSolver::SolveResult DiagonalLowRankSolver::Solve(
- const SolveData& data, std::vector<std::vector<DComplex>>& solutions,
+ const FullSolveData& data, std::vector<std::vector<DComplex>>& solutions,
double time, std::ostream* stat_stream) {
PrepareConstraints();
@@ -120,11 +120,12 @@ DiagonalLowRankSolver::SolveResult DiagonalLowRankSolver::Solve(
return result;
}
-void DiagonalLowRankSolver::CalculateNormPerDirection(const SolveData& data) {
+void DiagonalLowRankSolver::CalculateNormPerDirection(
+ const FullSolveData& data) {
std::vector<std::pair<float, size_t>> norm_sum_per_direction(NDirections());
for (size_t channel_block_index = 0;
channel_block_index != data.NChannelBlocks(); ++channel_block_index) {
- const SolveData::ChannelBlockData& cb_data =
+ const FullSolveData::ChannelBlockData& cb_data =
data.ChannelBlock(channel_block_index);
const size_t n_visibilities = cb_data.NVisibilities();
for (size_t direction = 0; direction != cb_data.NDirections();
@@ -147,7 +148,7 @@ void DiagonalLowRankSolver::CalculateNormPerDirection(const SolveData& data) {
}
double DiagonalLowRankSolver::ChiSquared(
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
std::vector<aocommon::MC2x2F>& v_residual, size_t direction,
const SolutionSpan& solutions) const {
using DComplex = std::complex<double>;
@@ -182,7 +183,7 @@ double DiagonalLowRankSolver::ChiSquared(
}
void DiagonalLowRankSolver::PerformIteration(
- size_t ch_block, const SolveData::ChannelBlockData& cb_data,
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
std::vector<MC2x2F>& v_residual, const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions, size_t iteration) {
// Fill v_residual
@@ -245,7 +246,7 @@ void AddToCorrelation(std::complex<float>& correlation_element,
}
void DiagonalLowRankSolver::SolveDirectionSolution(
- size_t ch_block, const SolveData::ChannelBlockData& cb_data,
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
const std::vector<aocommon::MC2x2F>& v_residual, size_t direction_index,
size_t solution_index, const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions) {
diff --git a/ddecal/gain_solvers/DiagonalLowRankSolver.h b/ddecal/gain_solvers/DiagonalLowRankSolver.h
index 9e9e98d5ada8074d8eb1f9bc01e3da5272207282..bf58625d836fcf62eb503124da83ccd47cb1e1bc 100644
--- a/ddecal/gain_solvers/DiagonalLowRankSolver.h
+++ b/ddecal/gain_solvers/DiagonalLowRankSolver.h
@@ -56,7 +56,7 @@ class DiagonalLowRankSolver final : public SolverBase {
n_low_rank_approximation_iterations),
n_power_iterations_(n_power_iterations) {}
- SolveResult Solve(const SolveData& data,
+ SolveResult Solve(const FullSolveData& data,
std::vector<std::vector<DComplex>>& solutions, double time,
std::ostream* stat_stream) override;
@@ -79,7 +79,7 @@ class DiagonalLowRankSolver final : public SolverBase {
private:
void PerformIteration(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
std::vector<aocommon::MC2x2F>& v_residual,
const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions, size_t iteration);
@@ -88,18 +88,18 @@ class DiagonalLowRankSolver final : public SolverBase {
* Calculates the chi-squared value, i.e. the weighted squared sum of the
* difference between the data and the corrected residual.
*/
- double ChiSquared(const SolveData::ChannelBlockData& cb_data,
+ double ChiSquared(const FullSolveData::ChannelBlockData& cb_data,
std::vector<aocommon::MC2x2F>& v_residual, size_t direction,
const SolutionSpan& solutions) const;
void SolveDirectionSolution(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
const std::vector<aocommon::MC2x2F>& v_residual,
size_t direction_index, size_t solution_index,
const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions);
- void CalculateNormPerDirection(const SolveData& data);
+ void CalculateNormPerDirection(const FullSolveData& data);
std::vector<size_t> direction_ordering_;
size_t n_low_rank_approximation_iterations_ = 25;
diff --git a/ddecal/gain_solvers/DiagonalSolver.cc b/ddecal/gain_solvers/DiagonalSolver.cc
index 1510c25766ea6d20f34187efc563c468739427c6..c31a557aae79d2641720dc0e5d6c5ac6102cfce1 100644
--- a/ddecal/gain_solvers/DiagonalSolver.cc
+++ b/ddecal/gain_solvers/DiagonalSolver.cc
@@ -20,7 +20,7 @@ namespace dp3 {
namespace ddecal {
DiagonalSolver::SolveResult DiagonalSolver::Solve(
- const SolveData& data, std::vector<std::vector<DComplex>>& solutions,
+ const FullSolveData& data, std::vector<std::vector<DComplex>>& solutions,
double time, std::ostream* stat_stream) {
assert(solutions.size() == NChannelBlocks());
@@ -65,7 +65,7 @@ DiagonalSolver::SolveResult DiagonalSolver::Solve(
[&](size_t start_block, size_t end_block, size_t thread_index) {
for (size_t ch_block = start_block; ch_block != end_block;
++ch_block) {
- const SolveData::ChannelBlockData& channel_block =
+ const FullSolveData::ChannelBlockData& channel_block =
data.ChannelBlock(ch_block);
std::vector<Matrix>& g_times_cs = thread_g_times_cs[thread_index];
@@ -112,7 +112,7 @@ DiagonalSolver::SolveResult DiagonalSolver::Solve(
}
void DiagonalSolver::PerformIteration(
- size_t ch_block, const SolveData::ChannelBlockData& cb_data,
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
std::vector<Matrix>& g_times_cs, std::vector<std::vector<Complex>>& vs,
const std::vector<DComplex>& solutions, SolutionTensor& next_solutions) {
for (size_t ant_and_pol = 0; ant_and_pol != NAntennas() * 2; ++ant_and_pol) {
@@ -213,7 +213,7 @@ static void Reset(std::vector<SolverBase::Complex>& vector, size_t size) {
}
void DiagonalSolver::InitializeModelMatrix(
- const SolveData::ChannelBlockData& channel_block_data,
+ const FullSolveData::ChannelBlockData& channel_block_data,
std::vector<Matrix>& g_times_cs,
std::vector<std::vector<Complex>>& vs) const {
assert(g_times_cs.empty() == vs.empty());
diff --git a/ddecal/gain_solvers/DiagonalSolver.h b/ddecal/gain_solvers/DiagonalSolver.h
index cf20b3ff3d7a2fa970a790bfba3391eda0238c3a..df4a38099575485095664ed7e32a29b0fd381a54 100644
--- a/ddecal/gain_solvers/DiagonalSolver.h
+++ b/ddecal/gain_solvers/DiagonalSolver.h
@@ -14,7 +14,7 @@ class DiagonalSolver final : public SolverBase {
public:
DiagonalSolver() : SolverBase() {}
- SolveResult Solve(const SolveData& data,
+ SolveResult Solve(const FullSolveData& data,
std::vector<std::vector<DComplex>>& solutions, double time,
std::ostream* stat_stream) override;
@@ -22,7 +22,7 @@ class DiagonalSolver final : public SolverBase {
private:
void PerformIteration(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
std::vector<Matrix>& g_times_cs,
std::vector<std::vector<Complex>>& vs,
const std::vector<DComplex>& solutions,
@@ -35,7 +35,7 @@ class DiagonalSolver final : public SolverBase {
* antenna visibilities in the corresponding channel block.
*/
void InitializeModelMatrix(
- const SolveData::ChannelBlockData& channel_block_data,
+ const FullSolveData::ChannelBlockData& channel_block_data,
std::vector<Matrix>& g_times_cs,
std::vector<std::vector<Complex>>& vs) const;
};
diff --git a/ddecal/gain_solvers/FullJonesSolver.cc b/ddecal/gain_solvers/FullJonesSolver.cc
index b101ffbed8e7b288245268361a7b6be740351342..2076ee7234340ab6e13b10aae8e738353229f007 100644
--- a/ddecal/gain_solvers/FullJonesSolver.cc
+++ b/ddecal/gain_solvers/FullJonesSolver.cc
@@ -16,7 +16,7 @@ namespace dp3 {
namespace ddecal {
FullJonesSolver::SolveResult FullJonesSolver::Solve(
- const SolveData& data, std::vector<std::vector<DComplex>>& solutions,
+ const FullSolveData& data, std::vector<std::vector<DComplex>>& solutions,
double time, std::ostream* stat_stream) {
// This algorithm is basically the same as the scalar algorithm,
// but visibility values are extended to 2x2 matrices and concatenated
@@ -89,7 +89,7 @@ FullJonesSolver::SolveResult FullJonesSolver::Solve(
[&](size_t start_block, size_t end_block, size_t thread_index) {
for (size_t ch_block = start_block; ch_block != end_block;
++ch_block) {
- const SolveData::ChannelBlockData& channel_block =
+ const FullSolveData::ChannelBlockData& channel_block =
data.ChannelBlock(ch_block);
std::vector<Matrix>& g_times_cs = thread_g_times_cs[thread_index];
@@ -137,7 +137,7 @@ FullJonesSolver::SolveResult FullJonesSolver::Solve(
}
void FullJonesSolver::PerformIteration(
- size_t ch_block, const SolveData::ChannelBlockData& cb_data,
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
std::vector<Matrix>& g_times_cs, std::vector<Matrix>& vs,
const std::vector<DComplex>& solutions, SolutionTensor& next_solutions) {
for (size_t ant = 0; ant != NAntennas(); ++ant) {
@@ -255,7 +255,7 @@ void FullJonesSolver::PerformIteration(
}
void FullJonesSolver::InitializeModelMatrix(
- const SolveData::ChannelBlockData& channel_block_data,
+ const FullSolveData::ChannelBlockData& channel_block_data,
std::vector<Matrix>& g_times_cs, std::vector<Matrix>& vs) const {
assert(g_times_cs.empty() == vs.empty());
if (g_times_cs.empty()) {
diff --git a/ddecal/gain_solvers/FullJonesSolver.h b/ddecal/gain_solvers/FullJonesSolver.h
index 66b279d1a1d853719a24a917272a6329bcecfa04..715f53053509eefc702dca875711b84fb59d53c4 100644
--- a/ddecal/gain_solvers/FullJonesSolver.h
+++ b/ddecal/gain_solvers/FullJonesSolver.h
@@ -12,7 +12,7 @@ namespace ddecal {
class FullJonesSolver final : public SolverBase {
public:
- SolveResult Solve(const SolveData& data,
+ SolveResult Solve(const FullSolveData& data,
std::vector<std::vector<DComplex>>& solutions, double time,
std::ostream* stat_stream) override;
@@ -20,7 +20,7 @@ class FullJonesSolver final : public SolverBase {
private:
void PerformIteration(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
std::vector<Matrix>& g_times_cs,
std::vector<Matrix>& vs,
const std::vector<DComplex>& solutions,
@@ -33,7 +33,7 @@ class FullJonesSolver final : public SolverBase {
* antenna visibilities in the corresponding channel block.
*/
void InitializeModelMatrix(
- const SolveData::ChannelBlockData& channel_block_data,
+ const FullSolveData::ChannelBlockData& channel_block_data,
std::vector<Matrix>& g_times_cs, std::vector<Matrix>& vs) const;
};
diff --git a/ddecal/gain_solvers/HybridSolver.cc b/ddecal/gain_solvers/HybridSolver.cc
index 57f87dfd2bc7ce6d0e32bf2e7913e0db2cd35d33..e0c41fd883f0adedd34570da847025da262c4015 100644
--- a/ddecal/gain_solvers/HybridSolver.cc
+++ b/ddecal/gain_solvers/HybridSolver.cc
@@ -22,8 +22,9 @@ void HybridSolver::AddSolver(std::unique_ptr<SolverBase> solver) {
}
HybridSolver::SolveResult HybridSolver::Solve(
- const SolveData& solve_data, std::vector<std::vector<DComplex>>& solutions,
- double time, std::ostream* stat_stream) {
+ const FullSolveData& solve_data,
+ std::vector<std::vector<DComplex>>& solutions, double time,
+ std::ostream* stat_stream) {
assert(!solvers_.empty());
size_t available_iterations = SolverBase::GetMaxIterations();
SolveResult result;
@@ -40,7 +41,8 @@ HybridSolver::SolveResult HybridSolver::Solve(
}
bool HybridSolver::RunSolver(SolverBase& solver, size_t& available_iterations,
- SolveResult& result, const SolveData& solve_data,
+ SolveResult& result,
+ const FullSolveData& solve_data,
std::vector<std::vector<DComplex>>& solutions,
double time, std::ostream* stat_stream) {
if (solver.GetMaxIterations() > available_iterations)
diff --git a/ddecal/gain_solvers/HybridSolver.h b/ddecal/gain_solvers/HybridSolver.h
index a7ad6d991403a92c50394099bd1c1dd861d29879..451f28e7089a24fec41ec2aee994dec31155a609 100644
--- a/ddecal/gain_solvers/HybridSolver.h
+++ b/ddecal/gain_solvers/HybridSolver.h
@@ -72,13 +72,13 @@ class HybridSolver final : public SolverBase {
}
/** @} */
- SolveResult Solve(const SolveData& solve_data,
+ SolveResult Solve(const FullSolveData& solve_data,
std::vector<std::vector<DComplex>>& solutions, double time,
std::ostream* stat_stream) override;
private:
bool RunSolver(SolverBase& solver, size_t& available_iterations,
- SolveResult& result, const SolveData& solve_data,
+ SolveResult& result, const FullSolveData& solve_data,
std::vector<std::vector<DComplex>>& solutions, double time,
std::ostream* stat_stream);
diff --git a/ddecal/gain_solvers/IterativeDiagonalSolver.cc b/ddecal/gain_solvers/IterativeDiagonalSolver.cc
index 82d0154b31516b7165a879d536ba22c603f4c801..b24b83894eb93e448b6ca0ef63cb2cc55bbeb641 100644
--- a/ddecal/gain_solvers/IterativeDiagonalSolver.cc
+++ b/ddecal/gain_solvers/IterativeDiagonalSolver.cc
@@ -16,7 +16,7 @@ namespace dp3 {
namespace ddecal {
IterativeDiagonalSolver::SolveResult IterativeDiagonalSolver::Solve(
- const SolveData& data, std::vector<std::vector<DComplex>>& solutions,
+ const FullSolveData& data, std::vector<std::vector<DComplex>>& solutions,
double time, std::ostream* stat_stream) {
PrepareConstraints();
@@ -84,7 +84,7 @@ IterativeDiagonalSolver::SolveResult IterativeDiagonalSolver::Solve(
}
void IterativeDiagonalSolver::PerformIteration(
- size_t ch_block, const SolveData::ChannelBlockData& cb_data,
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
std::vector<MC2x2F>& v_residual, const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions) {
// Fill v_residual
@@ -111,7 +111,7 @@ void IterativeDiagonalSolver::PerformIteration(
}
void IterativeDiagonalSolver::SolveDirection(
- size_t ch_block, const SolveData::ChannelBlockData& cb_data,
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
const std::vector<MC2x2F>& v_residual, size_t direction,
const std::vector<DComplex>& solutions, SolutionTensor& next_solutions) {
// Calculate this equation, given ant a:
diff --git a/ddecal/gain_solvers/IterativeDiagonalSolver.h b/ddecal/gain_solvers/IterativeDiagonalSolver.h
index 475254c4bb45d26cfca6ebb30cf5a69113c106b0..9c388843a75d8fddef6c0a2d83ea5fbed401d082 100644
--- a/ddecal/gain_solvers/IterativeDiagonalSolver.h
+++ b/ddecal/gain_solvers/IterativeDiagonalSolver.h
@@ -12,7 +12,7 @@ namespace ddecal {
class IterativeDiagonalSolver final : public SolverBase {
public:
- SolveResult Solve(const SolveData& data,
+ SolveResult Solve(const FullSolveData& data,
std::vector<std::vector<DComplex>>& solutions, double time,
std::ostream* stat_stream) override;
@@ -22,13 +22,13 @@ class IterativeDiagonalSolver final : public SolverBase {
private:
void PerformIteration(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
std::vector<aocommon::MC2x2F>& v_residual,
const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions);
void SolveDirection(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
const std::vector<aocommon::MC2x2F>& v_residual,
size_t direction, const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions);
diff --git a/ddecal/gain_solvers/IterativeFullJonesSolver.cc b/ddecal/gain_solvers/IterativeFullJonesSolver.cc
index f9996cd8ab6fc9e250e80770d8a19e53f2455a27..b0d951f54094ee9cdaacfa12680a1374bd5d42d3 100644
--- a/ddecal/gain_solvers/IterativeFullJonesSolver.cc
+++ b/ddecal/gain_solvers/IterativeFullJonesSolver.cc
@@ -15,7 +15,7 @@ namespace dp3 {
namespace ddecal {
IterativeFullJonesSolver::SolveResult IterativeFullJonesSolver::Solve(
- const SolveData& data, std::vector<std::vector<DComplex>>& solutions,
+ const FullSolveData& data, std::vector<std::vector<DComplex>>& solutions,
double time, std::ostream* stat_stream) {
PrepareConstraints();
@@ -83,7 +83,7 @@ IterativeFullJonesSolver::SolveResult IterativeFullJonesSolver::Solve(
}
void IterativeFullJonesSolver::PerformIteration(
- size_t ch_block, const SolveData::ChannelBlockData& cb_data,
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
std::vector<MC2x2F>& v_residual, const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions) {
// Fill v_residual
@@ -108,7 +108,7 @@ void IterativeFullJonesSolver::PerformIteration(
}
void IterativeFullJonesSolver::SolveDirection(
- size_t ch_block, const SolveData::ChannelBlockData& cb_data,
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
const std::vector<MC2x2F>& v_residual, size_t direction,
const std::vector<DComplex>& solutions, SolutionTensor& next_solutions) {
// Calculate this equation, given ant a:
@@ -180,8 +180,9 @@ void IterativeFullJonesSolver::SolveDirection(
template <bool Add>
void IterativeFullJonesSolver::AddOrSubtractDirection(
- const SolveData::ChannelBlockData& cb_data, std::vector<MC2x2F>& v_residual,
- size_t direction, const std::vector<DComplex>& solutions) {
+ const FullSolveData::ChannelBlockData& cb_data,
+ std::vector<MC2x2F>& v_residual, size_t direction,
+ const std::vector<DComplex>& solutions) {
constexpr size_t n_solution_polarizations = 4;
const size_t n_visibilities = cb_data.NVisibilities();
for (size_t vis_index = 0; vis_index != n_visibilities; ++vis_index) {
diff --git a/ddecal/gain_solvers/IterativeFullJonesSolver.h b/ddecal/gain_solvers/IterativeFullJonesSolver.h
index 98b6a84bb33702c122cd929a72f545dcde5b40a9..efb078e22996cd62da32b5ac418396cf86e5147d 100644
--- a/ddecal/gain_solvers/IterativeFullJonesSolver.h
+++ b/ddecal/gain_solvers/IterativeFullJonesSolver.h
@@ -16,7 +16,7 @@ namespace ddecal {
*/
class IterativeFullJonesSolver final : public SolverBase {
public:
- SolveResult Solve(const SolveData& data,
+ SolveResult Solve(const FullSolveData& data,
std::vector<std::vector<DComplex>>& solutions, double time,
std::ostream* stat_stream) override;
@@ -26,19 +26,19 @@ class IterativeFullJonesSolver final : public SolverBase {
private:
void PerformIteration(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
std::vector<aocommon::MC2x2F>& v_residual,
const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions);
template <bool Add>
- void AddOrSubtractDirection(const SolveData::ChannelBlockData& cb_data,
+ void AddOrSubtractDirection(const FullSolveData::ChannelBlockData& cb_data,
std::vector<aocommon::MC2x2F>& v_residual,
size_t direction,
const std::vector<DComplex>& solutions);
void SolveDirection(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
const std::vector<aocommon::MC2x2F>& v_residual,
size_t direction, const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions);
diff --git a/ddecal/gain_solvers/IterativeScalarSolver.cc b/ddecal/gain_solvers/IterativeScalarSolver.cc
index dfc0477877bb16160ed0d9ad6a82342f79f6b752..332f65e1ffbe918eac6bb9f9392c5b38bfae4e9e 100644
--- a/ddecal/gain_solvers/IterativeScalarSolver.cc
+++ b/ddecal/gain_solvers/IterativeScalarSolver.cc
@@ -16,7 +16,7 @@ namespace dp3 {
namespace ddecal {
IterativeScalarSolver::SolveResult IterativeScalarSolver::Solve(
- const SolveData& data, std::vector<std::vector<DComplex>>& solutions,
+ const FullSolveData& data, std::vector<std::vector<DComplex>>& solutions,
double time, std::ostream* stat_stream) {
PrepareConstraints();
@@ -84,7 +84,7 @@ IterativeScalarSolver::SolveResult IterativeScalarSolver::Solve(
}
void IterativeScalarSolver::PerformIteration(
- size_t ch_block, const SolveData::ChannelBlockData& cb_data,
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
std::vector<MC2x2F>& v_residual, const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions) {
// Fill v_residual
@@ -109,7 +109,7 @@ void IterativeScalarSolver::PerformIteration(
}
void IterativeScalarSolver::SolveDirection(
- size_t ch_block, const SolveData::ChannelBlockData& cb_data,
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
const std::vector<MC2x2F>& v_residual, size_t direction,
const std::vector<DComplex>& solutions, SolutionTensor& next_solutions) {
// Calculate this equation, given ant a:
@@ -169,8 +169,9 @@ void IterativeScalarSolver::SolveDirection(
template <bool Add>
void IterativeScalarSolver::AddOrSubtractDirection(
- const SolveData::ChannelBlockData& cb_data, std::vector<MC2x2F>& v_residual,
- size_t direction, const std::vector<DComplex>& solutions) {
+ const FullSolveData::ChannelBlockData& cb_data,
+ std::vector<MC2x2F>& v_residual, size_t direction,
+ const std::vector<DComplex>& solutions) {
const size_t n_visibilities = cb_data.NVisibilities();
for (size_t vis_index = 0; vis_index != n_visibilities; ++vis_index) {
const uint32_t antenna_1 = cb_data.Antenna1Index(vis_index);
diff --git a/ddecal/gain_solvers/IterativeScalarSolver.h b/ddecal/gain_solvers/IterativeScalarSolver.h
index 85156879943ffe8cbe9e9608026a0fc697d143ca..4765bb9177e259f8f6a37e79fd014e27d236211b 100644
--- a/ddecal/gain_solvers/IterativeScalarSolver.h
+++ b/ddecal/gain_solvers/IterativeScalarSolver.h
@@ -15,7 +15,7 @@ namespace ddecal {
*/
class IterativeScalarSolver final : public SolverBase {
public:
- SolveResult Solve(const SolveData& data,
+ SolveResult Solve(const FullSolveData& data,
std::vector<std::vector<DComplex>>& solutions, double time,
std::ostream* stat_stream) override;
@@ -25,19 +25,19 @@ class IterativeScalarSolver final : public SolverBase {
private:
void PerformIteration(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
std::vector<aocommon::MC2x2F>& v_residual,
const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions);
template <bool Add>
- void AddOrSubtractDirection(const SolveData::ChannelBlockData& cb_data,
+ void AddOrSubtractDirection(const FullSolveData::ChannelBlockData& cb_data,
std::vector<aocommon::MC2x2F>& v_residual,
size_t direction,
const std::vector<DComplex>& solutions);
void SolveDirection(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
const std::vector<aocommon::MC2x2F>& v_residual,
size_t direction, const std::vector<DComplex>& solutions,
SolutionTensor& next_solutions);
diff --git a/ddecal/gain_solvers/LBFGSSolver.cc b/ddecal/gain_solvers/LBFGSSolver.cc
index 5abe0cca717201f3bc8ce735e315cba3b172ed71..2cf1bbbe98e8831ffa4db3195e114054f5604cc2 100644
--- a/ddecal/gain_solvers/LBFGSSolver.cc
+++ b/ddecal/gain_solvers/LBFGSSolver.cc
@@ -34,13 +34,13 @@ namespace ddecal {
namespace {
struct lbfgs_fulljones_data {
- const SolveData::ChannelBlockData& cb_data;
+ const FullSolveData::ChannelBlockData& cb_data;
std::size_t n_antennas;
std::size_t n_directions;
double robust_nu;
std::size_t start_baseline;
std::size_t end_baseline;
- lbfgs_fulljones_data(const SolveData::ChannelBlockData& cb_data_,
+ lbfgs_fulljones_data(const FullSolveData::ChannelBlockData& cb_data_,
std::size_t n_antennas_, std::size_t n_directions_,
double robust_nu_)
: cb_data(cb_data_),
@@ -384,7 +384,7 @@ void ScalarGradient(double* unknowns, double* gradient, int n_unknowns,
}
void PerformIterationFull(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
const std::vector<SolverBase::DComplex>& solutions,
SolutionTensor& next_solutions,
std::size_t n_antennas, std::size_t n_directions,
@@ -420,7 +420,7 @@ void PerformIterationFull(size_t ch_block,
}
void PerformIterationDiagonal(
- size_t ch_block, const SolveData::ChannelBlockData& cb_data,
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
const std::vector<SolverBase::DComplex>& solutions,
SolutionTensor& next_solutions, std::size_t n_antennas,
std::size_t n_directions, double robust_nu, std::size_t max_iter,
@@ -454,7 +454,7 @@ void PerformIterationDiagonal(
}
void PerformIterationScalar(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
const std::vector<SolverBase::DComplex>& solutions,
SolutionTensor& next_solutions,
std::size_t n_antennas, std::size_t n_directions,
@@ -521,7 +521,7 @@ void LBFGSSolver::MergeSolutions(SolutionTensor& next_solutions,
}
LBFGSSolver::SolveResult LBFGSSolver::Solve(
- const SolveData& data, std::vector<std::vector<DComplex>>& solutions,
+ const FullSolveData& data, std::vector<std::vector<DComplex>>& solutions,
double time, std::ostream* stat_stream) {
assert(solutions.size() == NChannelBlocks());
@@ -537,7 +537,7 @@ LBFGSSolver::SolveResult LBFGSSolver::Solve(
{NChannelBlocks(), NAntennas(), NSolutions(), NSolutionPolarizations()});
for (size_t ch_block = 0; ch_block != NChannelBlocks(); ++ch_block) {
- const SolveData::ChannelBlockData& channel_block_data =
+ const FullSolveData::ChannelBlockData& channel_block_data =
data.ChannelBlock(ch_block);
// initialize with
// minibatches, size of parameters, size of data (baselines), history size,
diff --git a/ddecal/gain_solvers/LBFGSSolver.h b/ddecal/gain_solvers/LBFGSSolver.h
index 293ea92eeb492779884b4f0fd93f08386fcabd7e..aee4a21a89782dc33b695aac9176520e0943af59 100644
--- a/ddecal/gain_solvers/LBFGSSolver.h
+++ b/ddecal/gain_solvers/LBFGSSolver.h
@@ -47,7 +47,7 @@ class LBFGSSolver final : public SolverBase {
static void MergeSolutions(SolutionTensor& next_solutions, size_t ch_block,
const xt::xtensor<double, 1>& d_storage);
- SolveResult Solve(const SolveData& data,
+ SolveResult Solve(const FullSolveData& data,
std::vector<std::vector<DComplex>>& solutions, double time,
std::ostream* stat_stream) override;
diff --git a/ddecal/gain_solvers/ScalarSolver.cc b/ddecal/gain_solvers/ScalarSolver.cc
index 2eb1a6202aa80b0d0ba5d6aeeac21fb4292ba5b2..d06406369f2ec8f89e98f2edba2715ade4ebc74b 100644
--- a/ddecal/gain_solvers/ScalarSolver.cc
+++ b/ddecal/gain_solvers/ScalarSolver.cc
@@ -18,7 +18,7 @@ namespace dp3 {
namespace ddecal {
ScalarSolver::SolveResult ScalarSolver::Solve(
- const SolveData& data, std::vector<std::vector<DComplex>>& solutions,
+ const FullSolveData& data, std::vector<std::vector<DComplex>>& solutions,
double time, std::ostream* stat_stream) {
assert(solutions.size() == NChannelBlocks());
@@ -63,7 +63,7 @@ ScalarSolver::SolveResult ScalarSolver::Solve(
[&](size_t start_block, size_t end_block, size_t thread_index) {
for (size_t ch_block = start_block; ch_block != end_block;
++ch_block) {
- const SolveData::ChannelBlockData& channel_block =
+ const FullSolveData::ChannelBlockData& channel_block =
data.ChannelBlock(ch_block);
std::vector<Matrix>& g_times_cs = thread_g_times_cs[thread_index];
@@ -109,12 +109,10 @@ ScalarSolver::SolveResult ScalarSolver::Solve(
return result;
}
-void ScalarSolver::PerformIteration(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
- std::vector<Matrix>& g_times_cs,
- std::vector<Matrix>& vs,
- const std::vector<DComplex>& solutions,
- SolutionTensor& next_solutions) {
+void ScalarSolver::PerformIteration(
+ size_t ch_block, const FullSolveData::ChannelBlockData& cb_data,
+ std::vector<Matrix>& g_times_cs, std::vector<Matrix>& vs,
+ const std::vector<DComplex>& solutions, SolutionTensor& next_solutions) {
const size_t n_visibilities = cb_data.NVisibilities();
const size_t p1_to_p2[4] = {0, 2, 1, 3};
@@ -194,7 +192,7 @@ void ScalarSolver::PerformIteration(size_t ch_block,
}
void ScalarSolver::InitializeModelMatrix(
- const SolveData::ChannelBlockData& channel_block_data,
+ const FullSolveData::ChannelBlockData& channel_block_data,
std::vector<Matrix>& g_times_cs, std::vector<Matrix>& vs) const {
assert(g_times_cs.empty() == vs.empty());
if (g_times_cs.empty()) {
diff --git a/ddecal/gain_solvers/ScalarSolver.h b/ddecal/gain_solvers/ScalarSolver.h
index 3435ab5bc020c92b0be4bf74f11419c35d25dbea..de6aa2fb67ed6f894f5b73c68d04bb13f8d26901 100644
--- a/ddecal/gain_solvers/ScalarSolver.h
+++ b/ddecal/gain_solvers/ScalarSolver.h
@@ -12,7 +12,7 @@ namespace ddecal {
class ScalarSolver final : public SolverBase {
public:
- SolveResult Solve(const SolveData& data,
+ SolveResult Solve(const FullSolveData& data,
std::vector<std::vector<DComplex>>& solutions, double time,
std::ostream* stat_stream) override;
@@ -20,7 +20,7 @@ class ScalarSolver final : public SolverBase {
private:
void PerformIteration(size_t ch_block,
- const SolveData::ChannelBlockData& cb_data,
+ const FullSolveData::ChannelBlockData& cb_data,
std::vector<Matrix>& g_times_cs,
std::vector<Matrix>& vs,
const std::vector<DComplex>& solutions,
@@ -33,7 +33,7 @@ class ScalarSolver final : public SolverBase {
* antenna visibilities in the corresponding channel block.
*/
void InitializeModelMatrix(
- const SolveData::ChannelBlockData& channel_block_data,
+ const FullSolveData::ChannelBlockData& channel_block_data,
std::vector<Matrix>& g_times_cs, std::vector<Matrix>& vs) const;
};
diff --git a/ddecal/gain_solvers/SolveData.cc b/ddecal/gain_solvers/SolveData.cc
index 36f03cb2fb2a0cba9ee7a3d3d1373d000e586848..2732453dc28022f542ad2c6a88f2658d401f9876 100644
--- a/ddecal/gain_solvers/SolveData.cc
+++ b/ddecal/gain_solvers/SolveData.cc
@@ -14,12 +14,23 @@ using dp3::base::BDABuffer;
namespace dp3 {
namespace ddecal {
-SolveData::SolveData(const std::vector<base::DPBuffer>& buffers,
- const std::vector<std::string>& direction_names,
- size_t n_channel_blocks, size_t n_antennas,
- const std::vector<size_t>& n_solutions_per_direction,
- const std::vector<int>& antennas1,
- const std::vector<int>& antennas2)
+template <typename MatrixType>
+MatrixType ToSpecificMatrix(const std::complex<float>* data);
+template <>
+aocommon::MC2x2F ToSpecificMatrix(const std::complex<float>* data) {
+ return aocommon::MC2x2F(data);
+}
+template <>
+aocommon::MC2x2FDiag ToSpecificMatrix(const std::complex<float>* data) {
+ return aocommon::MC2x2FDiag(data[0], data[3]);
+}
+
+template <typename MatrixType>
+SolveData<MatrixType>::SolveData(
+ const std::vector<base::DPBuffer>& buffers,
+ const std::vector<std::string>& direction_names, size_t n_channel_blocks,
+ size_t n_antennas, const std::vector<size_t>& n_solutions_per_direction,
+ const std::vector<int>& antennas1, const std::vector<int>& antennas2)
: channel_blocks_(n_channel_blocks) {
std::vector<size_t> channel_begin(n_channel_blocks + 1, 0);
@@ -93,8 +104,8 @@ SolveData::SolveData(const std::vector<base::DPBuffer>& buffers,
const size_t channel_block_size = end_channel - first_channel;
for (size_t i = 0; i < channel_block_size; ++i) {
- cb_data.data_[vis_index + i] =
- aocommon::MC2x2F(&data(baseline, first_channel + i, 0));
+ cb_data.data_[vis_index + i] = ToSpecificMatrix<MatrixType>(
+ &data(baseline, first_channel + i, 0));
cb_data.antenna_indices_[vis_index + i] =
std::pair<uint32_t, uint32_t>(antenna1, antenna2);
}
@@ -118,7 +129,8 @@ SolveData::SolveData(const std::vector<base::DPBuffer>& buffers,
for (size_t i = 0; i < channel_block_size; ++i) {
cb_data.model_data_(direction, vis_index + i) =
- aocommon::MC2x2F(&model_data(baseline, first_channel + i, 0));
+ ToSpecificMatrix<MatrixType>(
+ &model_data(baseline, first_channel + i, 0));
cb_data.solution_map_(direction, vis_index + i) = solution_index;
}
@@ -133,10 +145,13 @@ SolveData::SolveData(const std::vector<base::DPBuffer>& buffers,
CountAntennaVisibilities(n_antennas);
}
-SolveData::SolveData(const BdaSolverBuffer& buffer, size_t n_channel_blocks,
- size_t n_directions, size_t n_antennas,
- const std::vector<int>& antennas1,
- const std::vector<int>& antennas2, bool with_weights)
+template <typename MatrixType>
+SolveData<MatrixType>::SolveData(const BdaSolverBuffer& buffer,
+ size_t n_channel_blocks, size_t n_directions,
+ size_t n_antennas,
+ const std::vector<int>& antennas1,
+ const std::vector<int>& antennas2,
+ bool with_weights)
: channel_blocks_(n_channel_blocks) {
// Count nr of visibilities
std::vector<size_t> counts(n_channel_blocks, 0);
@@ -184,8 +199,8 @@ SolveData::SolveData(const BdaSolverBuffer& buffer, size_t n_channel_blocks,
const std::complex<float>* data_ptr =
data_row.data + channel_start * data_row.n_correlations;
for (size_t i = 0; i != channel_block_size; ++i) {
- cb_data.data_[vis_index + i] =
- aocommon::MC2x2F(&data_ptr[i * data_row.n_correlations]);
+ cb_data.data_[vis_index + i] = ToSpecificMatrix<MatrixType>(
+ &data_ptr[i * data_row.n_correlations]);
cb_data.antenna_indices_[vis_index + i] =
std::pair<uint32_t, uint32_t>(antenna1, antenna2);
}
@@ -205,8 +220,9 @@ SolveData::SolveData(const BdaSolverBuffer& buffer, size_t n_channel_blocks,
model_data_row.data +
channel_start * model_data_row.n_correlations;
for (size_t i = 0; i != channel_block_size; ++i) {
- cb_data.model_data_(dir, vis_index + i) = aocommon::MC2x2F(
- &model_data_ptr[i * model_data_row.n_correlations]);
+ cb_data.model_data_(dir, vis_index + i) =
+ ToSpecificMatrix<MatrixType>(
+ &model_data_ptr[i * model_data_row.n_correlations]);
}
}
@@ -220,7 +236,8 @@ SolveData::SolveData(const BdaSolverBuffer& buffer, size_t n_channel_blocks,
cb_data.InitializeSolutionIndices(); // TODO replace!
}
-void SolveData::CountAntennaVisibilities(size_t n_antennas) {
+template <typename MatrixType>
+void SolveData<MatrixType>::CountAntennaVisibilities(size_t n_antennas) {
for (ChannelBlockData& cb_data : channel_blocks_) {
cb_data.antenna_visibility_counts_.assign(n_antennas, 0);
for (const std::pair<uint32_t, uint32_t>& a : cb_data.antenna_indices_) {
@@ -230,7 +247,8 @@ void SolveData::CountAntennaVisibilities(size_t n_antennas) {
}
}
-void SolveData::ChannelBlockData::InitializeSolutionIndices() {
+template <typename MatrixType>
+void SolveData<MatrixType>::ChannelBlockData::InitializeSolutionIndices() {
// This initializes the solution indices for direction-independent intervals
// TODO support DD intervals
n_solutions_.assign(NDirections(), 1);
@@ -239,5 +257,8 @@ void SolveData::ChannelBlockData::InitializeSolutionIndices() {
}
}
+template class SolveData<aocommon::MC2x2F>;
+template class SolveData<aocommon::MC2x2FDiag>;
+
} // namespace ddecal
} // namespace dp3
diff --git a/ddecal/gain_solvers/SolveData.h b/ddecal/gain_solvers/SolveData.h
index 88864029f4dedb0d80b77d4f0c80d3d602ffa831..6e790703253f7254e614b970d01215802c3dc39a 100644
--- a/ddecal/gain_solvers/SolveData.h
+++ b/ddecal/gain_solvers/SolveData.h
@@ -20,13 +20,15 @@ class BdaSolverBuffer;
/**
* Contains exactly the data required for solving: (weighted) data, model_data
* and the associated antennas for each visibility. In this class, the term
- * visibility refers to a 2x2 matrix containing the 4 polarizations.
+ * visibility refers to either a 2x2 diagonal or a full 2x2 matrix, containing 2
+ * or 4 polarizations respectively.
*/
+template <typename MatrixType = aocommon::MC2x2F>
class SolveData {
public:
class ChannelBlockData {
public:
- using const_iterator = std::vector<aocommon::MC2x2F>::const_iterator;
+ using const_iterator = typename std::vector<MatrixType>::const_iterator;
void Resize(size_t n_visibilities, size_t n_directions) {
data_.resize(n_visibilities);
@@ -70,11 +72,8 @@ class SolveData {
const uint32_t* SolutionMapData() const { return solution_map_.data(); }
const float& Weight(size_t index) const { return weights_(index, 0); }
- const aocommon::MC2x2F& Visibility(size_t index) const {
- return data_[index];
- }
- const aocommon::MC2x2F& ModelVisibility(size_t direction,
- size_t index) const {
+ const MatrixType& Visibility(size_t index) const { return data_[index]; }
+ const MatrixType& ModelVisibility(size_t direction, size_t index) const {
return model_data_(direction, index);
}
@@ -89,19 +88,19 @@ class SolveData {
}
private:
- friend class SolveData;
+ friend class SolveData<MatrixType>;
/**
* Initialize n_solutions_ and solution_map_.
*/
void InitializeSolutionIndices();
- std::vector<aocommon::MC2x2F> data_;
+ std::vector<MatrixType> data_;
// weights_(i, pol) contains the weight for data_[i][pol]. The vector will
// be left empty when the algorithm does not need the weights.
xt::xtensor<float, 2> weights_;
// model_data_(d, i) is the model data for direction d, element i
- xt::xtensor<aocommon::MC2x2F, 2> model_data_;
+ xt::xtensor<MatrixType, 2> model_data_;
// Element i contains the first and second antenna corresponding with
// data_[i] and model_data_(d, i)
std::vector<std::pair<uint32_t, uint32_t>> antenna_indices_;
@@ -162,14 +161,20 @@ class SolveData {
std::vector<ChannelBlockData> channel_blocks_;
};
-template <bool Add>
+/// Stores all 4 polarizations of the data.
+using FullSolveData = SolveData<aocommon::MC2x2F>;
+/// Stores only the 2 diagonal values of the data (e.g. XX/YY). Because the word
+/// "diagonal" is extensively used for the solution type, the term "Duo" is used
+/// for this.
+using DuoSolveData = SolveData<aocommon::MC2x2FDiag>;
+
+template <bool Add, typename MatrixType>
void DiagonalAddOrSubtractDirection(
- const SolveData::ChannelBlockData& cb_data,
- std::vector<aocommon::MC2x2F>& v_residual, size_t direction,
- size_t n_solutions, const std::vector<std::complex<double>>& solutions) {
+ const typename SolveData<MatrixType>::ChannelBlockData& cb_data,
+ std::vector<MatrixType>& v_residual, size_t direction, size_t n_solutions,
+ const std::vector<std::complex<double>>& solutions) {
using DComplex = std::complex<double>;
using Complex = std::complex<float>;
- using aocommon::MC2x2F;
const size_t n_visibilities = cb_data.NVisibilities();
for (size_t vis_index = 0; vis_index != n_visibilities; ++vis_index) {
const uint32_t antenna_1 = cb_data.Antenna1Index(vis_index);
@@ -179,17 +184,16 @@ void DiagonalAddOrSubtractDirection(
&solutions[(antenna_1 * n_solutions + solution_index) * 2];
const DComplex* solution_2 =
&solutions[(antenna_2 * n_solutions + solution_index) * 2];
- const Complex solution_1_0(solution_1[0]);
- const Complex solution_1_1(solution_1[1]);
- const Complex solution_2_0_conj(std::conj(solution_2[0]));
- const Complex solution_2_1_conj(std::conj(solution_2[1]));
-
- MC2x2F& data = v_residual[vis_index];
- const MC2x2F& model = cb_data.ModelVisibility(direction, vis_index);
- const MC2x2F contribution(solution_1_0 * model[0] * solution_2_0_conj,
- solution_1_0 * model[1] * solution_2_1_conj,
- solution_1_1 * model[2] * solution_2_0_conj,
- solution_1_1 * model[3] * solution_2_1_conj);
+
+ MatrixType& data = v_residual[vis_index];
+ const MatrixType& model = cb_data.ModelVisibility(direction, vis_index);
+ // Convert first to single precision to make calculation easier
+ const aocommon::MC2x2FDiag solution_a(static_cast<Complex>(solution_1[0]),
+ static_cast<Complex>(solution_1[1]));
+ const aocommon::MC2x2FDiag solution_b(static_cast<Complex>(solution_2[0]),
+ static_cast<Complex>(solution_2[1]));
+ const MatrixType contribution(solution_a * model *
+ HermTranspose(solution_b));
if constexpr (Add)
data += contribution;
else
@@ -197,6 +201,9 @@ void DiagonalAddOrSubtractDirection(
}
}
+extern template class SolveData<aocommon::MC2x2F>;
+extern template class SolveData<aocommon::MC2x2FDiag>;
+
} // namespace ddecal
} // namespace dp3
diff --git a/ddecal/gain_solvers/SolverBase.h b/ddecal/gain_solvers/SolverBase.h
index 6a780bd0f624b893b21f6237752b631d55a2e052..3c9910122df1ca9615facdabe3bbc273e57c74d3 100644
--- a/ddecal/gain_solvers/SolverBase.h
+++ b/ddecal/gain_solvers/SolverBase.h
@@ -13,12 +13,11 @@
#include "../constraints/Constraint.h"
#include "../linear_solvers/LLSSolver.h"
+#include "SolveData.h"
namespace dp3 {
namespace ddecal {
-class SolveData;
-
class SolverBase {
public:
typedef std::complex<double> DComplex;
@@ -201,7 +200,7 @@ class SolverBase {
* pol solutions.
* @param statStream Optional pointer to a stream for displaying statistics.
*/
- virtual SolveResult Solve(const SolveData& data,
+ virtual SolveResult Solve(const FullSolveData& data,
std::vector<std::vector<DComplex>>& solutions,
double time, std::ostream* statStream) = 0;
diff --git a/ddecal/test/unit/tSolveData.cc b/ddecal/test/unit/tSolveData.cc
index 25b4d3b4b212b499103dce4b08a1d83816473bbe..b1557fa6eae081f882d86930cc99c6e75e78d41a 100644
--- a/ddecal/test/unit/tSolveData.cc
+++ b/ddecal/test/unit/tSolveData.cc
@@ -16,7 +16,7 @@
using dp3::base::BDABuffer;
using dp3::base::DPBuffer;
-using ChannelBlockData = dp3::ddecal::SolveData::ChannelBlockData;
+using ChannelBlockData = dp3::ddecal::FullSolveData::ChannelBlockData;
namespace {
diff --git a/external/aocommon b/external/aocommon
index 2363df2242671d1835ca67952c16e3bde4a12148..0fe1970ad441f704225b6a64545e6c1c80053026 160000
--- a/external/aocommon
+++ b/external/aocommon
@@ -1 +1 @@
-Subproject commit 2363df2242671d1835ca67952c16e3bde4a12148
+Subproject commit 0fe1970ad441f704225b6a64545e6c1c80053026