DDECal.h

Go to the documentation of this file.

// DDECal.h: DP3 step class to calibrate direction dependent gains
// Copyright (C) 2023 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later
 
 
#ifndef DP3_STEPS_DDECAL_H_
#define DP3_STEPS_DDECAL_H_
 
#include <fstream>
#include <string>
#include <vector>
 
#include <aocommon/recursivefor.h>
 
#include <schaapcommon/h5parm/h5cache.h>
#include <schaapcommon/h5parm/jonesparameters.h>
 
#include "common/ParameterSet.h"
 
#include "ddecal/Settings.h"
#include "ddecal/SolutionWriter.h"
#include "ddecal/constraints/Constraint.h"
#include "ddecal/gain_solvers/SolverBase.h"
 
#include "MultiResultStep.h"
#include "ResultStep.h"
#include "UVWFlagger.h"
 
namespace dp3 {
namespace steps {
 
class DDECal : public Step {
 public:
  DDECal(const common::ParameterSet& parameterSet, const std::string& prefix);
 
  common::Fields getRequiredFields() const override {
    return kDataField | kFlagsField | kWeightsField | kUvwField;
  }
 
  common::Fields getProvidedFields() const override {
    return (settings_.subtract ||
            (settings_.only_predict && !settings_.keep_model_data))
               ? kDataField
               : common::Fields();
  }
 
  bool process(std::unique_ptr<base::DPBuffer> buffer) override;
 
  void checkMinimumVisibilities(size_t bufferIndex);
 
  void flagChannelBlock(size_t cbIndex, size_t bufferIndex);
 
  void doSolve();
 
  void finish() override;
 
  void updateInfo(const base::DPInfo&) override;
 
  void show(std::ostream&) const override;
 
  void showTimings(std::ostream&, double duration) const override;
 
 private:
  void initializeColumnReaders(const common::ParameterSet&,
                               const std::string& prefix);
  void initializeModelReuse();
  void initializeInitialSolutionsH5Parm(const common::ParameterSet& parset,
                                        const std::string& prefix);
  void initializeIDG(const common::ParameterSet& parset,
                     const std::string& prefix);
  void initializePredictSteps(const common::ParameterSet& parset,
                              const std::string& prefix);
 
  void setModelNextSteps(Step&, const std::string& direction,
                         const common::ParameterSet& parset,
                         const std::string& prefix) const;
 
  void doPrepare();
 
  void InitializeSolutions(size_t buffer_index);
 
  void WriteSolutions();
 
  void SumModels(size_t buffer_index);
 
  void ApplySolution(
      base::DPBuffer& buffer, size_t baseline, size_t channel,
      const std::vector<std::complex<double>>& channel_block_solutions) const;
 
  void CorrectAndSubtractModels(size_t buffer_index);
 
  xt::xtensor<std::complex<float>, 3> ReadJonesMatrixFromH5Parm(
      const base::Direction& direction, double timestamp,
      schaapcommon::h5parm::GainType gain_type,
      schaapcommon::h5parm::SolTab* first_soltab,
      schaapcommon::h5parm::SolTab* second_soltab);
 
  ddecal::Settings settings_;
 
  std::vector<std::vector<std::unique_ptr<base::DPBuffer>>> input_buffers_;
  xt::xtensor<bool, 5> original_flags_;
 
  double average_time_;
 
  std::vector<std::vector<std::vector<casacore::DComplex>>> solutions_;
  std::vector<size_t> n_iterations_;  // Number of iterations taken
  std::vector<size_t> n_approximating_iterations_;
 
  std::vector<std::vector<std::vector<ddecal::ConstraintResult>>>
      constraint_solutions_;
 
  std::unique_ptr<ddecal::SolutionWriter> solution_writer_;
 
  size_t requested_solution_interval_;
  size_t n_solution_intervals_;  
  size_t first_solution_index_;
  size_t n_channels_;
  std::vector<std::pair<size_t, size_t>> visibilities_in_interval_;
  std::vector<size_t> channel_block_start_;
  std::vector<double> channel_block_frequencies_;
  std::vector<std::vector<std::string>> patches_per_direction_;
  std::vector<std::string> direction_names_;
  std::vector<std::string> reused_direction_names_;
  std::vector<base::Direction> source_directions_;
 
  std::vector<int> antennas1_;
  std::vector<int> antennas2_;
  std::vector<double> weights_per_antenna_;
 
  UVWFlagger uvw_flag_step_;
  std::shared_ptr<ResultStep> data_result_step_;
  std::vector<std::shared_ptr<ModelDataStep>> steps_;
  std::vector<common::Fields> required_fields_;
  std::vector<std::shared_ptr<MultiResultStep>> result_steps_;
 
  bool store_solution_in_buffer_;
 
  std::unique_ptr<schaapcommon::h5parm::H5Parm> initial_solutions_;
  std::string initial_solutions_h5_parm_name;
  std::vector<std::string> initial_solutions_table_;
  std::vector<schaapcommon::h5parm::SolTab> solution_tables_;
  bool initial_solutions_are_full_jones_;
  schaapcommon::h5parm::JonesParameters::InterpolationType interpolation_type_;
  schaapcommon::h5parm::JonesParameters::MissingAntennaBehavior
      missing_antenna_behavior_;
  std::vector<schaapcommon::h5parm::GainType> gain_types_;
 
  common::NSTimer timer_;
  common::NSTimer predict_timer_;
  common::NSTimer solve_timer_;
  common::NSTimer write_timer_;
  std::mutex measures_mutex_;
  std::unique_ptr<ddecal::SolverBase> solver_;
  std::unique_ptr<std::ofstream> statistics_stream_;
};
 
}  // namespace steps
}  // namespace dp3
 
#endif