casacore.h

Go to the documentation of this file.

#ifndef PURIFY_CASACORE_H
#define PURIFY_CASACORE_H
 
#include "purify/config.h"
#include "purify/types.h"
#include <exception>
#include <map>
#include <memory>
#include <string>
#include <type_traits>
#include <casacore/ms/MeasurementSets/MeasurementSet.h>
#include <casacore/tables/TaQL/TableParse.h>
#include <casacore/tables/Tables/ArrayColumn.h>
#include <casacore/tables/Tables/ScalarColumn.h>
#include <casacore/tables/Tables/Table.h>
#include "purify/uvw_utilities.h"
#include <sopt/utilities.h>
 
namespace purify {
namespace casa {
 
template <class T>
Matrix<T> table_column(::casacore::Table const &table, std::string const &column,
                       std::string const &filter = "");
 
class MeasurementSet {
 public:
  class const_iterator;
  class ChannelWrapper;
  static std::string const default_filter;
  typedef Eigen::Array<t_real, 2, 1> Direction;
 
  MeasurementSet(std::string const filename)
      : filename_(filename),
        tables_(std::make_shared<std::map<std::string, ::casacore::Table>>()){};
  MeasurementSet(MeasurementSet const &c) : filename_(c.filename_), tables_(c.tables_){};
  std::string const &filename() const { return filename_; }
  MeasurementSet &filename(std::string const &filename);
 
  ::casacore::Table const &table(std::string const &name = "") const;
 
  template <class T>
  ::casacore::ScalarColumn<T> scalar_column(std::string const &col,
                                            std::string const &tabname = "") const {
    return get_column<T, ::casacore::ScalarColumn>(col, tabname);
  }
  template <class T>
  ::casacore::ArrayColumn<T> array_column(std::string const &col,
                                          std::string const &tabname = "") const {
    return get_column<T, ::casacore::ArrayColumn>(col, tabname);
  }
  void clear() { tables_ = std::make_shared<std::map<std::string, ::casacore::Table>>(); }
 
  template <class T>
  Matrix<T> column(std::string const &column, std::string const &filter = "") const {
    return table_column<T>(this->table(), column, filter);
  }
 
  std::size_t size() const;
 
  const_iterator begin(std::string const &filter = "") const;
  const_iterator end(std::string const &filter = "") const;
  ChannelWrapper operator[](t_uint i) const;
  ChannelWrapper operator[](std::tuple<t_uint, std::string> const &i) const;
  Direction direction(t_real tolerance = 1e-8, std::string const &filter = "") const;
  Direction::Scalar right_ascension(t_real tolerance = 1e-8, std::string const &filter = "") const {
    return direction(tolerance, filter)(0);
  }
  Direction::Scalar declination(t_real tolerance = 1e-8, std::string const &filter = "") const {
    return direction(tolerance, filter)(1);
  }
 
 private:
  template <class T, template <class> class TYPE>
  TYPE<T> get_column(std::string const &col, std::string const &tabname) const {
    return {table(tabname), col};
  }
 
  std::string filename_;
 
  std::shared_ptr<std::map<std::string, ::casacore::Table>> tables_;
};
 
namespace details {
template <class C>
Matrix<C> get_taql_array(::casacore::TaQLResult const &taql) {
  auto const col = ::casacore::ArrayColumn<C>(taql.table(), "R");
  auto const data_column = col.getColumn();
  auto const shape = col.shape(0);
  auto nsize =
      std::accumulate(shape.begin(), shape.end(), 1, [](ssize_t a, ssize_t b) { return a * b; });
  typedef Eigen::Matrix<C, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Matrix;
  return Matrix::Map(data_column.data(), col.nrow(), nsize);
}
template <class C>
Matrix<C> get_taql_scalar(::casacore::TaQLResult const &taql) {
  auto const col = ::casacore::ScalarColumn<C>(taql.table(), "R");
  auto const data_column = col.getColumn();
  typedef Eigen::Matrix<C, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Matrix;
  return Matrix::Map(data_column.data(), col.nrow(), 1);
}
 
template <class T>
Matrix<T> table_column(::casacore::Table const &table, std::string const &column,
                       std::string const &filter, std::integral_constant<bool, true> const &) {
  if (table.nrow() == 0) return Matrix<T>::Zero(0, 1);
  auto const taql_table = ::casacore::tableCommand(
      "USING STYLE PYTHON SELECT " + column + " as R FROM $1 " + filter, table);
  auto const vtable = taql_table.table();
  if (vtable.nrow() == 0) return Matrix<T>(0, 1);
 
  switch (vtable.tableDesc().columnDesc("R").trueDataType()) {
#define PURIFY_MACRO(NAME, TYPE) \
  case ::casacore::Tp##NAME:     \
    return details::get_taql_scalar<::casacore::TYPE>(taql_table).template cast<T>();
    PURIFY_MACRO(Complex, Complex);
    PURIFY_MACRO(DComplex, DComplex);
#undef PURIFY_MACRO
#define PURIFY_MACRO(NAME, TYPE)  \
  case ::casacore::TpArray##NAME: \
    return details::get_taql_array<::casacore::TYPE>(taql_table).template cast<T>();
    PURIFY_MACRO(Complex, Complex);
    PURIFY_MACRO(DComplex, DComplex);
#undef PURIFY_MACRO
  default:
    break;
  }
 
  throw std::runtime_error("Array type is not handled");
  return Matrix<T>::Zero(0, 1);
}
 
template <class T>
Matrix<T> table_column(::casacore::Table const &table, std::string const &column,
                       std::string const &filter, std::integral_constant<bool, false> const &) {
  if (table.nrow() == 0) return Matrix<T>::Zero(0, 1);
  auto const taql_table = ::casacore::tableCommand(
      "USING STYLE PYTHON SELECT " + column + " as R FROM $1 " + filter, table);
  auto const vtable = taql_table.table();
  if (vtable.nrow() == 0) return Matrix<T>(0, 1);
 
  switch (vtable.tableDesc().columnDesc("R").trueDataType()) {
#define PURIFY_MACRO(NAME, TYPE) \
  case ::casacore::Tp##NAME:     \
    return details::get_taql_scalar<::casacore::TYPE>(taql_table).template cast<T>();
    PURIFY_MACRO(Bool, Bool);
    PURIFY_MACRO(Char, Char);
    PURIFY_MACRO(UChar, uChar);
    PURIFY_MACRO(Short, Short);
    PURIFY_MACRO(UShort, uShort);
    PURIFY_MACRO(Int, Int);
    PURIFY_MACRO(UInt, uInt);
    PURIFY_MACRO(Float, Float);
    PURIFY_MACRO(Double, Double);
#undef PURIFY_MACRO
#define PURIFY_MACRO(NAME, TYPE)  \
  case ::casacore::TpArray##NAME: \
    return details::get_taql_array<::casacore::TYPE>(taql_table).template cast<T>();
    PURIFY_MACRO(Bool, Bool);
    PURIFY_MACRO(Char, Char);
    PURIFY_MACRO(UChar, uChar);
    PURIFY_MACRO(Short, Short);
    PURIFY_MACRO(UShort, uShort);
    PURIFY_MACRO(Int, Int);
    PURIFY_MACRO(UInt, uInt);
    PURIFY_MACRO(Float, Float);
    PURIFY_MACRO(Double, Double);
#undef PURIFY_MACRO
  default:
    break;
  }
 
  throw std::runtime_error("Array type is not handled");
  return Matrix<T>::Zero(0, 1);
}
}  // namespace details
 
template <class T>
Matrix<T> table_column(::casacore::Table const &table, std::string const &column,
                       std::string const &filter) {
  return details::table_column<T>(table, column, filter,
                                  std::integral_constant<bool, sopt::is_complex<T>::value>());
}
 
class MeasurementSet::ChannelWrapper {
 public:
  enum class Sigma { OVERALL, SPECTRUM };
  ChannelWrapper(t_uint channel, MeasurementSet const &ms, std::string const &filter = "")
      : ms_(ms), filter_(filter), channel_(channel) {}
 
  t_uint channel() const { return channel_; }
#define PURIFY_MACRO(NAME, INDEX)                                                 \
                                              \
  Vector<t_real> raw_##NAME() const {                                             \
    return table_column<t_real>(ms_.table(), "UVW[" #INDEX "]", filter());        \
  }                                                                               \
                         \
  Vector<t_real> lambda_##NAME() const {                                          \
    return (raw_##NAME().array() * frequencies().array()).matrix() / constant::c; \
  }
  PURIFY_MACRO(u, 0);
  PURIFY_MACRO(v, 1);
  PURIFY_MACRO(w, 2);
#undef PURIFY_MACRO
  t_uint size() const;
 
  Vector<t_int> field_ids() const { return ms_.column<t_int>("FIELD_ID", filter()); }
  Vector<t_int> data_desc_id() const { return ms_.column<t_int>("DATA_DESC_ID", filter()); }
 
  Direction direction(t_real tolerance = 1e-8) const { return ms_.direction(tolerance, filter()); }
  Direction::Scalar right_ascension(t_real tolerance = 1e-8) const {
    return ms_.right_ascension(tolerance, filter());
  }
  Direction::Scalar declination(t_real tolerance = 1e-8) const {
    return ms_.declination(tolerance, filter());
  }
 
  Vector<t_real> raw_frequencies() const { return raw_spectral_window("CHAN_FREQ"); }
 
#define PURIFY_MACRO(NAME)                                                                     \
                                                                      \
  Vector<t_complex> NAME(std::string const &col = "DATA") const {                              \
    return ms_.column<t_complex>(stokes(#NAME, index(col)), filter());                         \
  }                                                                                            \
                                           \
  Vector<t_real> w##NAME(Sigma const &col = Sigma::OVERALL) const {                            \
    return table_column<t_real>(                                                               \
        ms_.table(), stokes(#NAME, col == Sigma::OVERALL ? "SIGMA" : index("SIGMA_SPECTRUM")), \
        filter());                                                                             \
  }
  // Stokes parameters
  PURIFY_MACRO(I);
  PURIFY_MACRO(Q);
  PURIFY_MACRO(U);
  PURIFY_MACRO(V);
  // Circular correlations
  PURIFY_MACRO(LL);
  PURIFY_MACRO(RR);
  PURIFY_MACRO(LR);
  PURIFY_MACRO(RL);
  // Linear correlations
  PURIFY_MACRO(XX);
  PURIFY_MACRO(YY);
  PURIFY_MACRO(XY);
  PURIFY_MACRO(YX);
 
#undef PURIFY_MACRO
 
  Vector<t_real> frequencies() const { return joined_spectral_window("CHAN_FREQ"); }
  Vector<t_real> width() const { return joined_spectral_window("CHAN_WIDTH"); }
  Vector<t_real> effective_noise_bandwidth() const {
    return joined_spectral_window("EFFECTIVE_BW");
  }
  Vector<t_real> resolution() const { return joined_spectral_window("RESOLUTION"); }
 
  bool is_valid() const;
 
 protected:
  std::string filter() const;
  std::string index(std::string const &variable = "") const;
  std::string stokes(std::string const &pol, std::string const &column = "DATA") const;
 
  Vector<t_real> raw_spectral_window(std::string const &column) const;
  Vector<t_real> joined_spectral_window(std::string const &column) const;
 
  MeasurementSet ms_;
  std::string const filter_;
  t_uint channel_;
};
 
class MeasurementSet::const_iterator {
 public:
  typedef MeasurementSet::ChannelWrapper value_type;
  typedef std::shared_ptr<value_type const> pointer;
  typedef value_type const &reference;
  typedef t_int difference_type;
 
  const_iterator(t_int channel, MeasurementSet const &ms, std::string const &filter = "")
      : channel_(channel),
        ms_(ms),
        filter_(filter),
        wrapper_(new value_type(channel_, ms_, filter_)){};
 
  pointer operator->() const { return wrapper_; }
  reference operator*() const { return *wrapper_; }
  const_iterator &operator++();
  const_iterator operator++(int);
  const_iterator &operator+=(difference_type n);
  const_iterator &operator-=(difference_type n) { return operator+=(-n); }
  const_iterator operator+(difference_type n) const {
    return const_iterator(channel_ + n, ms_, filter_);
  }
  const_iterator operator-(difference_type n) const {
    return const_iterator(channel_ - n, ms_, filter_);
  }
  bool operator>(const_iterator const &c) const;
  bool operator>=(const_iterator const &c) const;
  bool operator<(const_iterator const &c) const { return not operator>=(c); }
  bool operator<=(const_iterator const &c) const { return not operator>(c); }
  bool operator==(const_iterator const &c) const;
  bool operator!=(const_iterator const &c) const { return not operator==(c); }
 
  bool same_measurement_set(const_iterator const &c) const {
    return &ms_.table() == &c.ms_.table();
  }
 
 protected:
  difference_type channel_;
  MeasurementSet ms_;
  std::string const filter_;
  std::shared_ptr<value_type> wrapper_;
};
utilities::vis_params read_measurementset(std::string const &filename, const stokes pol = stokes::I,
                                          const std::vector<t_int> &channels = std::vector<t_int>(),
                                          std::string const &filter = "");
utilities::vis_params read_measurementset(MeasurementSet const &ms_file,
                                          const stokes pol = stokes::I,
                                          const std::vector<t_int> &channels = std::vector<t_int>(),
                                          std::string const &filter = "");
utilities::vis_params read_measurementset(std::string const &filename,
                                          const utilities::vis_params &uv1, const stokes pol,
                                          const std::vector<t_int> &channels,
                                          std::string const &filter);
utilities::vis_params read_measurementset(std::vector<std::string> const &filename,
                                          const stokes pol = stokes::I,
                                          const std::vector<t_int> &channels = std::vector<t_int>(),
                                          std::string const &filter = "");
 
t_real average_frequency(const purify::casa::MeasurementSet &ms_file, std::string const &filter,
                         const std::vector<t_int> &channels);
std::vector<utilities::vis_params> read_measurementset_channels(std::string const &filename,
                                                                const stokes pol = stokes::I,
                                                                const t_int &channel_width = 1,
                                                                std::string const &filter = "");
 
inline MeasurementSet::const_iterator operator+(MeasurementSet::const_iterator::difference_type n,
                                                MeasurementSet::const_iterator const &c) {
  return c.operator+(n);
}
inline MeasurementSet::const_iterator operator-(MeasurementSet::const_iterator::difference_type n,
                                                MeasurementSet::const_iterator const &c) {
  return c.operator-(n);
}
}  // namespace casa
}  // namespace purify
 
#endif