/usr/include/purify/casacore.h is in libpurify-dev 2.0.0-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 | #ifndef PURIFY_CASACORE_H
#define PURIFY_CASACORE_H
#include "purify/config.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 <sopt/utilities.h>
#include "purify/types.h"
#include "purify/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 = "");
//! Interface around measurement sets
class MeasurementSet {
public:
//! Iterates over channels
class const_iterator;
class ChannelWrapper;
//! Default filter specifying which data to accept
static std::string const default_filter;
//! Type for (RA, DEC) direction
typedef Eigen::Array<t_real, 2, 1> Direction;
//! Constructs the interface around a given measurement set
MeasurementSet(std::string const filename)
: filename_(filename),
tables_(std::make_shared<std::map<std::string, ::casacore::Table>>()){};
//! Shallow measurement set copy
MeasurementSet(MeasurementSet const &c) : filename_(c.filename_), tables_(c.tables_){};
//! Filename of the measurement set
std::string const &filename() const { return filename_; }
//! Set new filename
MeasurementSet &filename(std::string const &filename);
//! \brief Gets table or subtable
::casacore::Table const &table(std::string const &name = "") const;
//! Gets scalar column from table
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);
}
//! \brief Gets array column from table
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);
}
//! Clear memory
void clear() { tables_ = std::make_shared<std::map<std::string, ::casacore::Table>>(); }
//! Data from a column
template <class T>
Matrix<T> column(std::string const &column, std::string const &filter = "") const {
return table_column<T>(this->table(), column, filter);
}
//! Number of channels in the measurement set
std::size_t size() const;
//! Iterates over channels
const_iterator begin(std::string const &filter = "") const;
//! Iterates over channels
const_iterator end(std::string const &filter = "") const;
//! Returns wrapper over specific channel
ChannelWrapper operator[](t_uint i) const;
//! Returns wrapper over specific channel
ChannelWrapper operator[](std::tuple<t_uint, std::string> const &i) const;
//! Direction (RA, DEC) in radians
Direction direction(t_real tolerance = 1e-8, std::string const &filter = "") const;
//! Right ascention in radians
Direction::Scalar right_ascension(t_real tolerance = 1e-8, std::string const &filter = "") const {
return direction(tolerance, filter)(0);
}
//! Declination in radians
Direction::Scalar declination(t_real tolerance = 1e-8, std::string const &filter = "") const {
return direction(tolerance, filter)(1);
}
private:
//! Gets stokes of given array/object
template <class T, template <class> class TYPE>
TYPE<T> get_column(std::string const &col, std::string const &tabname) const {
return {table(tabname), col};
}
//! Name associated with the measurement set
std::string filename_;
//! Holds tables that have been opened
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, "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, "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>(vtable).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>(vtable).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>(vtable).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>(vtable).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);
}
}
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:
//! Possible locations for SIGMA
enum class Sigma { OVERALL, SPECTRUM };
enum class polarization { I, Q, U, V, LL, RR, RL, LR, XX, YY, XY, YX };
ChannelWrapper(t_uint channel, MeasurementSet const &ms, std::string const &filter = "")
: ms_(ms), filter_(filter), channel_(channel) {}
//! Channel this object is associated with
t_uint channel() const { return channel_; }
#define PURIFY_MACRO(NAME, INDEX) \
/** Stokes component in meters **/ \
Vector<t_real> raw_##NAME() const { \
return table_column<t_real>(ms_.table(), "UVW[" #INDEX "]", filter()); \
} \
/** \brief U scaled to purify values of wavelengths **/ \
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
//! Number of rows in a channel
t_uint size() const;
//! FIELD_ID from table MAIN
Vector<t_int> field_ids() const { return ms_.column<t_int>("FIELD_ID", filter()); }
//! DATA_DESC_ID from table MAIN
Vector<t_int> data_desc_id() const { return ms_.column<t_int>("DATA_DESC_ID", filter()); }
//! Direction (RA, DEC) in radian
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());
}
//! Frequencies from SPECTRAL_WINDOW for a this channel
Vector<t_real> raw_frequencies() const { return raw_spectral_window("CHAN_FREQ"); }
#define PURIFY_MACRO(NAME) \
/** Stokes component */ \
Vector<t_complex> NAME(std::string const &col = "DATA") const { \
return ms_.column<t_complex>(stokes(#NAME, index(col)), filter()); \
} \
/** Standard deviation for the Stokes component */ \
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
//! Frequencies for each valid measurement
Vector<t_real> frequencies() const { return joined_spectral_window("CHAN_FREQ"); }
//! Channel width for each valid measurement
Vector<t_real> width() const { return joined_spectral_window("CHAN_WIDTH"); }
//! Effective noise band-width width for each valid measurement
Vector<t_real> effective_noise_bandwidth() const {
return joined_spectral_window("EFFECTIVE_BW");
}
//! Effective spectral resolution for each valid measurement
Vector<t_real> resolution() const { return joined_spectral_window("RESOLUTION"); }
//! Check if channel has any data
bool is_valid() const;
protected:
//! Composes filter for current channel
std::string filter() const;
//! Composes "<variable>[<channel>,]"
std::string index(std::string const &variable = "") const;
//! Computes given stokes polarization
std::string stokes(std::string const &pol, std::string const &column = "DATA") const;
//! column[<channel>] from SPECTRAL_WINDOW
Vector<t_real> raw_spectral_window(std::string const &column) const;
//! column[<channel>] from SPECTRAL_WINDOW joined to each measurement
Vector<t_real> joined_spectral_window(std::string const &column) const;
//! Table over which to operate
MeasurementSet ms_;
//! Filter to apply to searches
std::string const filter_;
//! Channel index
t_uint channel_;
};
class MeasurementSet::const_iterator {
public:
//! Convenience wrapper to access values associated with a single channel
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); }
//! True if iterating over the same measurement set
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_;
};
//! Read measurement set into vis_params structure
utilities::vis_params read_measurementset(std::string const &filename,
const MeasurementSet::ChannelWrapper::polarization pol
= MeasurementSet::ChannelWrapper::polarization::I,
const std::vector<t_int> &channels = std::vector<t_int>(),
std::string const &filter = "");
//! Return average frequency over channels
t_real average_frequency(const purify::casa::MeasurementSet &ms_file, std::string const &filter,
const std::vector<t_int> &channels);
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);
}
}
}
#endif
|