DiagonalLowRankSolver.h

Go to the documentation of this file.

// Copyright (C) 2024 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later
 
#ifndef DDECAL_DIAGONAL_LOW_RANK_SOLVER_H_
#define DDECAL_DIAGONAL_LOW_RANK_SOLVER_H_
 
#include "SolverBase.h"
#include "SolveData.h"
 
namespace dp3::ddecal {
 
float DominantEigenPair(const xt::xtensor<std::complex<float>, 2>& matrix,
                        xt::xtensor<std::complex<float>, 1>& eigen_vector,
                        size_t n_iterations);
 
float DominantEigenPairNear(const xt::xtensor<std::complex<float>, 2>& matrix,
                            xt::xtensor<std::complex<float>, 1>& eigen_vector,
                            size_t n_iterations);
 
class DiagonalLowRankSolver final : public SolverBase {
 public:
  DiagonalLowRankSolver(size_t n_low_rank_approximation_iterations,
                        size_t n_power_iterations)
      : SolverBase(),
        n_low_rank_approximation_iterations_(
            n_low_rank_approximation_iterations),
        n_power_iterations_(n_power_iterations) {}
 
  SolveResult Solve(const FullSolveData& data,
                    std::vector<std::vector<DComplex>>& solutions,
                    double time) override;
 
  size_t NSolutionPolarizations() const override { return 2; }
 
  void SetNPowerIterations(size_t n_power_iterations) {
    n_power_iterations_ = n_power_iterations;
  }
 
  void SetNLowRankApproximationIterations(size_t n_iterations) {
    n_low_rank_approximation_iterations_ = n_iterations;
  }
 
 private:
  void PerformIteration(size_t ch_block,
                        const FullSolveData::ChannelBlockData& cb_data,
                        std::vector<aocommon::MC2x2F>& v_residual,
                        const std::vector<DComplex>& solutions,
                        SolutionTensor& next_solutions, size_t iteration);
 
  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 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 FullSolveData& data);
 
  std::vector<size_t> direction_ordering_;
  size_t n_low_rank_approximation_iterations_ = 25;
  size_t n_power_iterations_ = 10;
};
 
}  // namespace dp3::ddecal
 
#endif  // DDECAL_DIAGONAL_LOW_RANK_SOLVER_H_