Loading...
Searching...
No Matches
UvwFile.h
Go to the documentation of this file.
#ifndef CASACORE_UVW_FILE_H_
#define CASACORE_UVW_FILE_H_
#include "BufferedColumnarFile.h"
#include <algorithm>
#include <cassert>
#include <limits>
#include <ostream>
#include <stdexcept>
#include <string>
namespace casacore {
class UvwFile {
public:
UvwFile() noexcept = default;
UvwFile(UvwFile&& source) noexcept
: file_(std::move(source.file_)),
n_rows_(source.n_rows_),
rows_per_block_(source.rows_per_block_),
active_block_(source.active_block_),
reference_antenna_(source.reference_antenna_),
start_antenna_2_(source.start_antenna_2_),
n_antennas_(source.n_antennas_),
block_uvws_(std::move(source.block_uvws_)),
block_is_changed_(source.block_is_changed_) {
source.n_rows_ = 0;
source.rows_per_block_ = 0;
source.active_block_ = 0;
source.reference_antenna_ = 0;
source.start_antenna_2_ = 0;
source.n_antennas_ = 0;
source.block_uvws_.clear();
source.block_is_changed_ = false;
}
~UvwFile() noexcept { Close(); }
UvwFile& operator=(UvwFile&& rhs) {
Close();
file_ = std::move(rhs.file_);
n_rows_ = rhs.n_rows_;
rows_per_block_ = rhs.rows_per_block_;
active_block_ = rhs.active_block_;
reference_antenna_ = rhs.reference_antenna_;
start_antenna_2_ = rhs.start_antenna_2_;
n_antennas_ = rhs.n_antennas_;
block_uvws_ = std::move(rhs.block_uvws_);
block_is_changed_ = rhs.block_is_changed_;
rhs.n_rows_ = 0;
rhs.rows_per_block_ = 0;
rhs.active_block_ = 0;
rhs.reference_antenna_ = 0;
rhs.start_antenna_2_ = 0;
rhs.n_antennas_ = 0;
rhs.block_uvws_.clear();
rhs.block_is_changed_ = false;
return *this;
}
static UvwFile CreateNew(const std::string& filename) {
return UvwFile(filename);
}
static UvwFile OpenExisting(const std::string& filename) {
return UvwFile(filename, true);
}
void WriteUvw(uint64_t row, size_t antenna1, size_t antenna2,
const double* uvw) {
assert(file_.IsOpen());
if (row > n_rows_) {
throw std::runtime_error(
"Uvw data must be written in order (writing row " +
std::to_string(row) + ", after writing " + std::to_string(n_rows_) +
" rows)");
}
// The row/block is zero when there's not yet a full block written.
if (rows_per_block_ == 0) {
if (row == 0) {
reference_antenna_ = antenna1;
start_antenna_2_ = antenna2;
n_antennas_ = std::max(antenna1, antenna2) + 1;
block_uvws_.resize(n_antennas_, kUnsetPosition);
block_uvws_[reference_antenna_] = {0.0, 0.0, 0.0};
} else if (antenna1 == reference_antenna_ &&
antenna2 == start_antenna_2_) {
// This baseline is the first baseline of a new block, so the block size
// can be determined
rows_per_block_ = n_rows_;
n_antennas_ = block_uvws_.size();
WriteHeader();
ActivateBlock(1);
} else {
n_antennas_ = std::max({antenna1 + 1, antenna2 + 1, n_antennas_});
block_uvws_.resize(n_antennas_, kUnsetPosition);
}
} else {
const uint64_t block = row / rows_per_block_;
ActivateBlock(block);
}
if (antenna1 != antenna2) {
if (antenna1 == reference_antenna_) {
// baseline = a2 - a1 with a1 = 0
const std::array<double, 3> ant2_uvw{uvw[0], uvw[1], uvw[2]};
StoreOrCheck(antenna2, ant2_uvw);
} else if (antenna2 == reference_antenna_) {
// baseline = a2 - a1 with a2 = 0
const std::array<double, 3> ant1_uvw{-uvw[0], -uvw[1], -uvw[2]};
StoreOrCheck(antenna1, ant1_uvw);
} else if (IsSet(antenna1)) {
// baseline = a2 - a1. Given a1:
// a2 = baseline + a1
const std::array<double, 3> ant1_uvw = block_uvws_[antenna1];
const std::array<double, 3> ant2_uvw{
uvw[0] + ant1_uvw[0], uvw[1] + ant1_uvw[1], uvw[2] + ant1_uvw[2]};
StoreOrCheck(antenna2, ant2_uvw);
} else if (IsSet(antenna2)) {
// baseline = a2 - a1. Given a2:
// a1 = a2 - baseline
const std::array<double, 3> ant2_uvw = block_uvws_[antenna2];
const std::array<double, 3> ant1_uvw{
ant2_uvw[0] - uvw[0], ant2_uvw[1] - uvw[1], ant2_uvw[2] - uvw[2]};
StoreOrCheck(antenna1, ant1_uvw);
} else {
throw std::runtime_error(
"Baselines are written in a non-ordered way: they need to be "
"ordered either by antenna 1 or by antenna 2");
}
} else {
// In the special case of a file with only auto-correlations, the
// positions of the antennas can not be established in this case, it's
// also not necessary, because a zero uvw can be returned for
// auto-correlations. However, it will cause StoreOrCheck() to be
// never called, and therefore the block is never written to the file, and
// upon read the nr. of rows can not be determined. So, if this is an
// auto-correlation that has not yet been written, mark the block as
// changed so it gets written.
block_is_changed_ = block_is_changed_ || !IsSet(antenna1);
}
n_rows_ = std::max(n_rows_, row + 1);
}
void ReadUvw(uint64_t row, size_t antenna1, size_t antenna2, double* uvw) {
assert(file_.IsOpen());
if (row >= n_rows_ || antenna1 >= n_antennas_ || antenna2 >= n_antennas_) {
throw std::runtime_error(
"Invalid read for Uvw data: row " + std::to_string(row) +
", baseline (" + std::to_string(antenna1) + ", " +
std::to_string(antenna2) + ") was requested. File has only " +
std::to_string(n_rows_) + " rows with " +
std::to_string(n_antennas_) + " antennas.");
}
if (rows_per_block_ != 0) {
const uint64_t block = row / rows_per_block_;
ActivateBlock(block);
}
uvw[0] = block_uvws_[antenna2][0] - block_uvws_[antenna1][0];
uvw[1] = block_uvws_[antenna2][1] - block_uvws_[antenna1][1];
uvw[2] = block_uvws_[antenna2][2] - block_uvws_[antenna1][2];
}
void Close() {
if (file_.IsOpen()) {
// This handles two special cases: i) if only one block of visibilities
// is written, no repetition of baseline would have been identified yet.
// In that case, we now know the block size. ii) in case only a single
// auto-correlation (one antenna) is written without any
// cross-correlations, the compressed file will remain empty. We then have
// to identify how many rows are really in the MS, which is done by
// setting rows_per_block_ to the nr of rows. When reading such an MS, the
// fact that n_antennas_ == 1 is then a trigger to use rows_per_block_ as
// nrows, instead of the size of the compressed file.
if (rows_per_block_ == 0) {
rows_per_block_ = n_rows_;
n_antennas_ = block_uvws_.size();
WriteHeader();
} else if (n_antennas_ == 1) {
rows_per_block_ = n_rows_;
WriteHeader();
}
if (block_is_changed_) {
WriteActiveBlock();
}
file_.Close();
}
}
uint64_t NRows() const { return n_rows_; }
std::string Filename() const { return file_.Filename(); }
private:
UvwFile(const std::string& filename)
: file_(BufferedColumnarFile::CreateNew(filename, kHeaderSize,
sizeof(double) * 3)) {}
UvwFile(const std::string& filename, bool /*open existing*/)
: file_(BufferedColumnarFile::OpenExisting(filename, kHeaderSize)) {
ReadHeader();
active_block_ = std::numeric_limits<uint64_t>::max();
if (n_antennas_ > 1) {
if (file_.NRows() % (n_antennas_ - 1) != 0) {
throw std::runtime_error(
"Uvw file has an incorrect number of rows (" +
std::to_string(file_.NRows()) + ", expecting multiple of " +
std::to_string(n_antennas_ - 1) + "): file corrupted?");
}
const uint64_t n_blocks = file_.NRows() / (n_antennas_ - 1);
n_rows_ = n_blocks * rows_per_block_;
} else if (n_antennas_ == 1) {
n_rows_ = rows_per_block_;
}
if (n_rows_ > 0 && n_antennas_ <= reference_antenna_)
throw std::runtime_error(
"Invalid combination of values for n_antenna and reference antenna "
"in file: file damaged?");
// Setting the size of block_uvws_ here, saves an extra size check in
// ActivateBlock().
block_uvws_.assign(n_antennas_, kUnsetPosition);
}
void StoreOrCheck(size_t antenna, const std::array<double, 3>& antenna_uvw) {
if (IsSet(antenna)) {
if (!AreNear(block_uvws_[antenna], antenna_uvw)) {
std::ostringstream msg;
msg << "Inconsistent UVW value written for antenna " << antenna
<< ": old value is " << UvwAsString(block_uvws_[antenna])
<< ", new value is " << UvwAsString(antenna_uvw) << ".";
throw std::runtime_error(msg.str());
}
} else {
if (block_uvws_.size() <= antenna)
block_uvws_.resize(antenna + 1, kUnsetPosition);
block_uvws_[antenna] = antenna_uvw;
block_is_changed_ = true;
}
}
void ActivateBlock(size_t block) {
if (block != active_block_) {
if (block_is_changed_) {
WriteActiveBlock();
}
active_block_ = block;
ReadActiveBlock();
}
}
void ReadActiveBlock() {
const uint64_t block_start_row = (n_antennas_ - 1) * active_block_;
if (block_start_row < file_.NRows()) {
block_uvws_.clear();
for (size_t antenna = 0; antenna != n_antennas_; ++antenna) {
if (antenna != reference_antenna_) {
const uint64_t row = antenna < reference_antenna_
? block_start_row + antenna
: block_start_row + antenna - 1;
std::array<double, 3>& uvw = block_uvws_.emplace_back();
file_.Read(row, 0, uvw.data(), 3);
} else {
block_uvws_.emplace_back(std::array<double, 3>{0.0, 0.0, 0.0});
}
}
} else {
std::fill(block_uvws_.begin(), block_uvws_.end(), kUnsetPosition);
block_uvws_[reference_antenna_] = {0.0, 0.0, 0.0};
}
}
void WriteActiveBlock() {
if (block_uvws_.size() != n_antennas_)
throw std::runtime_error("Trying to write an incomplete UVW block");
const uint64_t block_start_row = (n_antennas_ - 1) * active_block_;
for (size_t antenna = 0; antenna != n_antennas_; ++antenna) {
if (antenna != reference_antenna_) {
const uint64_t row = antenna < reference_antenna_
? block_start_row + antenna
: block_start_row + antenna - 1;
file_.Write(row, 0, block_uvws_[antenna].data(), 3);
}
}
block_is_changed_ = false;
}
void ReadHeader() {
unsigned char data[kHeaderSize];
file_.ReadHeader(data);
if (!std::equal(data, data + 8, kMagicHeaderTag)) {
throw std::runtime_error(
"The UVW columnar file header not have the expected tag for UVW "
"columns: the measurement set may be damaged");
}
rows_per_block_ = reinterpret_cast<uint64_t&>(data[8]);
reference_antenna_ = reinterpret_cast<uint64_t&>(data[16]);
n_antennas_ = reinterpret_cast<uint64_t&>(data[24]);
}
void WriteHeader() {
unsigned char data[kHeaderSize];
std::copy_n(kMagicHeaderTag, 8, data);
reinterpret_cast<uint64_t&>(data[8]) = rows_per_block_;
reinterpret_cast<uint64_t&>(data[16]) = reference_antenna_;
reinterpret_cast<uint64_t&>(data[24]) = n_antennas_;
file_.WriteHeader(data);
}
bool IsSet(size_t antenna) const {
return block_uvws_.size() > antenna &&
block_uvws_[antenna] != kUnsetPosition;
}
static bool AreNear(std::array<double, 3> a, std::array<double, 3> b) {
return AreNear(a[0], b[0]) && AreNear(a[1], b[1]) && AreNear(a[2], b[2]);
}
static bool AreNear(double a, double b) {
const double magnitude = std::max({1e-5, std::fabs(a), std::fabs(b)});
return (std::fabs(a - b) / magnitude) < 1e-5;
}
static std::string UvwAsString(const std::array<double, 3>& uvw) {
std::ostringstream str;
str << "[" << uvw[0] << ", " << uvw[1] << ", " << uvw[2] << "]";
return str.str();
}
constexpr static size_t kHeaderSize = 32;
constexpr static const char kMagicHeaderTag[8] = "Uvw-col";
constexpr static std::array<double, 3> kUnsetPosition = {
std::numeric_limits<double>::max(), std::numeric_limits<double>::max(),
std::numeric_limits<double>::max()};
BufferedColumnarFile file_;
uint64_t n_rows_ = 0;
uint64_t rows_per_block_ = 0;
uint64_t active_block_ = 0;
size_t reference_antenna_ = 0;
// This value is used to determine the first baseline in the data, which is
// the baseline (reference_antenna_, start_antenna_2_).
size_t start_antenna_2_ = 0;
size_t n_antennas_ = 0;
// UVW for each antenna in the block
std::vector<std::array<double, 3>> block_uvws_;
bool block_is_changed_ = false;
};
} // namespace casacore
#endif
BufferedColumnarFile.h
casacore::UvwFile::Filename
std::string Filename() const
Definition UvwFile.h:242
casacore::UvwFile::kHeaderSize
static constexpr size_t kHeaderSize
The header: char[8] "Uvw-col\0" (=kMagicHeaderTag) uint64_t rows_per_block uint64_t reference_antenna...
Definition UvwFile.h:398
casacore::UvwFile::active_block_
uint64_t active_block_
Definition UvwFile.h:423
casacore::UvwFile::OpenExisting
static UvwFile OpenExisting(const std::string &filename)
Open an already existing UVW file from disk with the given filename.
Definition UvwFile.h:107
casacore::UvwFile::kMagicHeaderTag
static constexpr const char kMagicHeaderTag[8]
Definition UvwFile.h:399
casacore::UvwFile::n_antennas_
size_t n_antennas_
Definition UvwFile.h:428
casacore::UvwFile::CreateNew
static UvwFile CreateNew(const std::string &filename)
Create a new UVW file on disk with the given filename.
Definition UvwFile.h:100
casacore::UvwFile::rows_per_block_
uint64_t rows_per_block_
A "block" is a contiguous number of baselines that together form one timestep.
Definition UvwFile.h:422
casacore::UvwFile::WriteUvw
void WriteUvw(uint64_t row, size_t antenna1, size_t antenna2, const double *uvw)
Write a single row to the column.
Definition UvwFile.h:116
casacore::UvwFile::block_is_changed_
bool block_is_changed_
Definition UvwFile.h:431
casacore::UvwFile::Close
void Close()
Definition UvwFile.h:215
casacore::UvwFile::file_
BufferedColumnarFile file_
Definition UvwFile.h:403
casacore::UvwFile::n_rows_
uint64_t n_rows_
Number of rows in the Uvw column.
Definition UvwFile.h:411
casacore::UvwFile::AreNear
static bool AreNear(std::array< double, 3 > a, std::array< double, 3 > b)
Definition UvwFile.h:378
casacore::UvwFile::block_uvws_
std::vector< std::array< double, 3 > > block_uvws_
UVW for each antenna in the block.
Definition UvwFile.h:430
casacore::UvwFile::NRows
uint64_t NRows() const
Definition UvwFile.h:241
casacore::UvwFile::AreNear
static bool AreNear(double a, double b)
Definition UvwFile.h:381
casacore::UvwFile::StoreOrCheck
void StoreOrCheck(size_t antenna, const std::array< double, 3 > &antenna_uvw)
If this block does not have a value for the specified antenna, store the uvw value for it.
Definition UvwFile.h:288
casacore::UvwFile::IsSet
bool IsSet(size_t antenna) const
Definition UvwFile.h:374
casacore::UvwFile::ReadHeader
void ReadHeader()
Definition UvwFile.h:352
casacore::UvwFile::UvwFile
UvwFile(const std::string &filename)
Create a new file on disk.
Definition UvwFile.h:248
casacore::UvwFile::kUnsetPosition
static constexpr std::array< double, 3 > kUnsetPosition
Definition UvwFile.h:400
casacore::UvwFile::operator=
UvwFile & operator=(UvwFile &&rhs)
Definition UvwFile.h:75
casacore::UvwFile::reference_antenna_
size_t reference_antenna_
Definition UvwFile.h:424
casacore::UvwFile::UvwAsString
static std::string UvwAsString(const std::array< double, 3 > &uvw)
Definition UvwFile.h:385
casacore::UvwFile::ReadActiveBlock
void ReadActiveBlock()
Definition UvwFile.h:316
casacore::UvwFile::UvwFile
UvwFile(const std::string &filename, bool)
Open an existing file from disk.
Definition UvwFile.h:256
casacore::UvwFile::UvwFile
UvwFile() noexcept=default
casacore::UvwFile::start_antenna_2_
size_t start_antenna_2_
This value is used to determine the first baseline in the data, which is the baseline (reference_ante...
Definition UvwFile.h:427
casacore::UvwFile::ReadUvw
void ReadUvw(uint64_t row, size_t antenna1, size_t antenna2, double *uvw)
Read a single row.
Definition UvwFile.h:196
casacore::UvwFile::WriteActiveBlock
void WriteActiveBlock()
Definition UvwFile.h:337
casacore::UvwFile::~UvwFile
~UvwFile() noexcept
Definition UvwFile.h:73
casacore::UvwFile::WriteHeader
void WriteHeader()
Definition UvwFile.h:365
casacore::UvwFile::ActivateBlock
void ActivateBlock(size_t block)
Definition UvwFile.h:305
casacore
this file contains all the compiler specific defines
Definition mainpage.dox:28
casacore::move
void move(TYPE *target, int npixels) const
casacore::BufferedColumnarFile
VarBufferedColumnarFile< 100 *1024 > BufferedColumnarFile
Definition BufferedColumnarFile.h:340
std
Define real & complex conjugation for non-complex types and put comparisons into std namespace.
Definition Complex.h:350