From 99f7bf84bd2d549c3083f86b0966f978bf9deb5b Mon Sep 17 00:00:00 2001 From: Maik Nijhuis <maik.nijhuis@triopsys.nl> Date: Mon, 27 May 2024 12:15:23 +0000 Subject: [PATCH] AST-1540 Merge baseline selection code --- CMakeLists.txt | 1 - base/BaselineSelection.cc | 137 ++++++++++++++++++++++++++++++++----- base/BaselineSelection.h | 33 +++++++-- common/BaselineSelect.cc | 139 -------------------------------------- common/BaselineSelect.h | 77 --------------------- steps/MSBDAReader.cc | 1 - steps/MSReader.cc | 57 +++++----------- 7 files changed, 165 insertions(+), 280 deletions(-) delete mode 100644 common/BaselineSelect.cc delete mode 100644 common/BaselineSelect.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 870d0fe48..4d34b1603 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -378,7 +378,6 @@ add_library( # 'Common' files are included in all executables and libraries. add_library( Common OBJECT - common/BaselineSelect.cc common/ClusterDesc.cc common/DataConvert.cc common/Fields.cc diff --git a/base/BaselineSelection.cc b/base/BaselineSelection.cc index a37f663fb..c447bcee6 100644 --- a/base/BaselineSelection.cc +++ b/base/BaselineSelection.cc @@ -6,22 +6,34 @@ #include "BaselineSelection.h" +#include <vector> + #include <boost/algorithm/string.hpp> -#include <casacore/casa/Utilities/Regex.h> #include <casacore/casa/Arrays/Matrix.h> +#include <casacore/casa/Utilities/Regex.h> +#include <casacore/casa/version.h> +#include <casacore/ms/MeasurementSets/MeasurementSet.h> +#include <casacore/ms/MeasurementSets/MSAntenna.h> +#include <casacore/ms/MeasurementSets/MSAntennaColumns.h> +#include <casacore/ms/MSSel/MSAntennaGram.h> +#include <casacore/tables/Tables/ScaColDesc.h> +#include <casacore/tables/Tables/SetupNewTab.h> +#include <casacore/tables/Tables/Table.h> + +#include <aocommon/logger.h> -#include "../common/BaselineSelect.h" #include "../common/ParameterSet.h" #include "../common/ParameterValue.h" #include "../common/StreamUtil.h" -#include <vector> - -#include <aocommon/logger.h> - using casacore::IPosition; using casacore::Matrix; +using casacore::MS; +using casacore::MSAntenna; +using casacore::MSAntennaParse; +using casacore::SetupNewTable; +using casacore::Table; using dp3::common::operator<<; using aocommon::Logger; @@ -29,6 +41,37 @@ using aocommon::Logger; namespace dp3 { namespace base { +LogAntennaParseErrors::LogAntennaParseErrors() + : old_handler_(MSAntennaParse::thisMSAErrorHandler) { + // The new handler logs all errors as warnings and does not throw exceptions. + class ErrorHandler : public casacore::MSSelectionErrorHandler { + public: + ErrorHandler() = default; + ~ErrorHandler() override = default; + void reportError(const char* token, + const casacore::String message) override { + Logger::Warn << message << token << '\n'; + } + }; + + // This syntax works both when ErrorHandlerPointer is a raw pointer + // (casacore < 3.1.2) and when it's a smart pointer. + MSAntennaParse::thisMSAErrorHandler = ErrorHandlerPointer(new ErrorHandler()); +} + +LogAntennaParseErrors::~LogAntennaParseErrors() { +#if CASACORE_MAJOR_VERSION < 3 || \ + (CASACORE_MAJOR_VERSION == 3 && \ + (CASACORE_MINOR_VERSION == 0 || \ + (CASACORE_MINOR_VERSION == 1 && CASACORE_PATCH_VERSION < 2))) + // In casacore < 3.1.2 thisMSAErrorHandler is a raw pointer, + // From casacore 3.1.2. it's a CountedPtr or another smart pointer. + delete MSAntennaParse::thisMSAErrorHandler; +#endif + + MSAntennaParse::thisMSAErrorHandler = old_handler_; +} + BaselineSelection::BaselineSelection() {} BaselineSelection::BaselineSelection(const common::ParameterSet& parset, @@ -125,17 +168,79 @@ void BaselineSelection::handleBL(Matrix<bool>& selectBL, info.antennaNames())); } else { // Specified in casacore's MSSelection format. - std::ostringstream os; - const casacore::Matrix<bool> sel = common::BaselineSelect::convert( - info.antennaNames(), info.antennaPos(), info.getAnt1(), info.getAnt2(), - itsStrBL, os); - // Show possible messages about unknown stations. - if (!os.str().empty()) { - Logger::Warn << os.str(); - } - assert(sel.nrow() == selectBL.nrow()); - selectBL = selectBL && sel; + selectBL = selectBL && HandleMsSelection(info); + } +} + +Matrix<bool> BaselineSelection::HandleMsSelection(const DPInfo& info) const { + const casacore::String& antenna1_string = MS::columnName(MS::ANTENNA1); + const casacore::String& antenna2_string = MS::columnName(MS::ANTENNA2); + + // Create a temporary MSAntenna table in memory for parsing purposes. + SetupNewTable antenna_setup(casacore::String(), + MSAntenna::requiredTableDesc(), Table::New); + Table antenna_table(antenna_setup, Table::Memory, info.antennaNames().size()); + MSAntenna ms_antenna(antenna_table); + casacore::MSAntennaColumns ms_antenna_columns(ms_antenna); + for (size_t i = 0; i < info.antennaNames().size(); ++i) { + ms_antenna_columns.name().put(i, info.antennaNames()[i]); + ms_antenna_columns.positionMeas().put(i, info.antennaPos()[i]); + } + + // Create a temporary table holding the antenna numbers of the baselines. + casacore::TableDesc table_description; + table_description.addColumn(casacore::ScalarColumnDesc<int>(antenna1_string)); + table_description.addColumn(casacore::ScalarColumnDesc<int>(antenna2_string)); + SetupNewTable setup(casacore::String(), table_description, Table::New); + Table table(setup, Table::Memory, info.nbaselines()); + casacore::ScalarColumn<int> column1(table, antenna1_string); + casacore::ScalarColumn<int> column2(table, antenna2_string); + for (size_t i = 0; i < info.nbaselines(); ++i) { + column1.put(i, info.getAnt1()[i]); + column2.put(i, info.getAnt2()[i]); } + + // Do the selection using the temporary tables. + casacore::TableExprNode antenna_node1 = table.col(antenna1_string); + casacore::TableExprNode antenna_node2 = table.col(antenna2_string); + + casacore::Vector<int> selected_antennas1; + casacore::Vector<int> selected_antennas2; + { + // Overwrite the error handler to ignore errors for unknown antennas. + dp3::base::LogAntennaParseErrors ignore_antenna_errors; + + // Parse the selection. + // 'selected_antennas[12]' will contain the selected antenna indices. + // 'selected_baselines' becomes an n x 2 matrix with the antenna indices + // for each selected baseline. + Matrix<int> selected_baselines; + casacore::TableExprNode selection_node = + casacore::msAntennaGramParseCommand( + antenna_table, antenna_node1, antenna_node2, itsStrBL, + selected_antennas1, selected_antennas2, selected_baselines); + + // msAntennaGramParseCommand may put negative indices (with unknown + // semantics) into 'selected_antennas[12]' and 'selected_baselines'. + // -> Apply 'selection_node' and extract the correct antenna indices. + Table selection_table = table(selection_node); + selected_antennas1 = + casacore::ScalarColumn<int>(selection_table, antenna1_string) + .getColumn(); + selected_antennas2 = + casacore::ScalarColumn<int>(selection_table, antenna2_string) + .getColumn(); + } + + // Convert selected_antennas[12] to a selection matrix. + Matrix<bool> selection(info.nantenna(), info.nantenna(), false); + for (size_t bl = 0; bl < selected_antennas1.size(); ++bl) { + const int a1 = selected_antennas1[bl]; + const int a2 = selected_antennas2[bl]; + selection(a1, a2) = true; + selection(a2, a1) = true; + } + return selection; } Matrix<bool> BaselineSelection::handleBLVector( diff --git a/base/BaselineSelection.h b/base/BaselineSelection.h index 2c904a5f4..23562301a 100644 --- a/base/BaselineSelection.h +++ b/base/BaselineSelection.h @@ -1,18 +1,20 @@ // BaselineSelection.h: Class to handle the baseline selection -// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy) +// Copyright (C) 2024 ASTRON (Netherlands Institute for Radio Astronomy) // SPDX-License-Identifier: GPL-3.0-or-later /// @file /// @brief Class to handle the baseline selection /// @author Ger van Diepen -#ifndef DPPP_BASELINESELECTION_H -#define DPPP_BASELINESELECTION_H +#ifndef DP3_BASELINESELECTION_H_ +#define DP3_BASELINESELECTION_H_ #include <dp3/base/DPInfo.h> #include <casacore/casa/Arrays/Vector.h> #include <casacore/casa/Arrays/Matrix.h> +#include <casacore/ms/MSSel/MSAntennaParse.h> +#include <casacore/ms/MSSel/MSSelectionErrorHandler.h> namespace dp3 { namespace common { @@ -22,6 +24,24 @@ class ParameterValue; namespace base { +/// RAII object for temporarily overriding MSAntennaParse::thisMSAErrorHandler. +class LogAntennaParseErrors { + public: + /// Constructor. Set MSAntennaParse::thisMSAErrorHandler to a handler + /// which forwards all errors as warnings to the Logger. + LogAntennaParseErrors(); + + /// Destructor. Restores the original MSAntennaParse::thisMSAErrorHandler + /// and cleans up the temporary installed handler. + ~LogAntennaParseErrors(); + + private: + // Different casacore versions use different (smart) pointer types. + using ErrorHandlerPointer = + decltype(casacore::MSAntennaParse::thisMSAErrorHandler); + ErrorHandlerPointer old_handler_; +}; + /// \brief Class containing a few static functions to parse a baseline selection /// string. class BaselineSelection { @@ -62,6 +82,9 @@ class BaselineSelection { /// Convert the baseline selection string. void handleBL(casacore::Matrix<bool>& selectBL, const DPInfo& info) const; + /// Handle an MSSelection string. + casacore::Matrix<bool> HandleMsSelection(const DPInfo& info) const; + /// Handle a vector of baseline specifications. casacore::Matrix<bool> handleBLVector( const common::ParameterValue& pvBL, @@ -73,8 +96,8 @@ class BaselineSelection { /// Handle the baseline length selection. void handleLength(casacore::Matrix<bool>& selectBL, const DPInfo& info) const; - string itsStrBL; - string itsCorrType; + std::string itsStrBL; + std::string itsCorrType; std::vector<double> itsRangeBL; }; diff --git a/common/BaselineSelect.cc b/common/BaselineSelect.cc deleted file mode 100644 index 8857af7ab..000000000 --- a/common/BaselineSelect.cc +++ /dev/null @@ -1,139 +0,0 @@ -// BaselineSelect.cc: Convert MSSelection baseline string to a Matrix -// -// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy) -// SPDX-License-Identifier: GPL-3.0-or-later - -#include "BaselineSelect.h" - -#include <cassert> - -#include <casacore/ms/MeasurementSets/MeasurementSet.h> -#include <casacore/ms/MeasurementSets/MSAntenna.h> -#include <casacore/ms/MeasurementSets/MSAntennaColumns.h> -#include <casacore/ms/MSSel/MSSelection.h> -#include <casacore/ms/MSSel/MSAntennaParse.h> -#include <casacore/ms/MSSel/MSSelectionErrorHandler.h> -#include <casacore/ms/MSSel/MSAntennaGram.h> -#include <casacore/tables/Tables/Table.h> -#include <casacore/tables/Tables/SetupNewTab.h> -#include <casacore/tables/Tables/TableRecord.h> -#include <casacore/tables/TaQL/TableParse.h> -#include <casacore/tables/Tables/ScalarColumn.h> -#include <casacore/tables/Tables/ScaColDesc.h> -#include <casacore/measures/Measures/MPosition.h> -#include <casacore/casa/Arrays/Matrix.h> -#include <casacore/casa/Arrays/Vector.h> -#include <casacore/casa/version.h> - -#include <dp3/common/Types.h> - -using casacore::Matrix; -using casacore::MPosition; -using casacore::MSAntenna; -using casacore::MSAntennaColumns; -using casacore::MSAntennaParse; -using casacore::MSSelectionErrorHandler; -using casacore::ScalarColumn; -using casacore::ScalarColumnDesc; -using casacore::SetupNewTable; -using casacore::Table; -using casacore::TableDesc; -using casacore::TableExprNode; - -namespace dp3 { -namespace common { - -casacore::Matrix<bool> BaselineSelect::convert( - const std::vector<std::string>& names, const std::vector<MPosition>& pos, - const std::vector<int>& antenna1, const std::vector<int>& antenna2, - const std::string& baseline_selection, std::ostream& os) { - assert(names.size() == pos.size()); - assert(antenna1.size() == antenna2.size()); - - // Create a temporary MSAntenna table in memory for parsing purposes. - SetupNewTable antNew(casacore::String(), MSAntenna::requiredTableDesc(), - Table::New); - Table anttab(antNew, Table::Memory, names.size()); - MSAntenna msant(anttab); - MSAntennaColumns antcol(msant); - antcol.name().putColumn(casacore::Vector<casacore::String>(names)); - for (size_t i = 0; i < pos.size(); ++i) { - antcol.positionMeas().put(i, pos[i]); - } - // Create a temporary table holding the antenna numbers of the baselines. - TableDesc td; - td.addColumn(ScalarColumnDesc<int>("ANTENNA1")); - td.addColumn(ScalarColumnDesc<int>("ANTENNA2")); - SetupNewTable tabNew(casacore::String(), td, Table::New); - Table tab(tabNew, Table::Memory, antenna1.size()); - ScalarColumn<int> ac1(tab, "ANTENNA1"); - ScalarColumn<int> ac2(tab, "ANTENNA2"); - ac1.putColumn(casacore::Vector<casacore::Int>(antenna1)); - ac2.putColumn(casacore::Vector<casacore::Int>(antenna2)); - - // Do the selection using the temporary tables. - TableExprNode a1(tab.col("ANTENNA1")); - TableExprNode a2(tab.col("ANTENNA2")); - return convert(anttab, a1, a2, baseline_selection, os); -} - -Matrix<bool> BaselineSelect::convert(Table& anttab, TableExprNode& a1Node, - TableExprNode& a2Node, - const string& baselineSelection, - std::ostream& os) { - // Overwrite the error handler to ignore errors for unknown antennas. - // First construct MSSelection, because it resets the error handler. - casacore::Vector<int> selectedAnts1; - casacore::Vector<int> selectedAnts2; - Matrix<int> selectedBaselines; - auto curHandler = MSAntennaParse::thisMSAErrorHandler; -#if CASACORE_MAJOR_VERSION < 3 || \ - (CASACORE_MAJOR_VERSION == 3 && \ - (CASACORE_MINOR_VERSION == 0 || \ - (CASACORE_MINOR_VERSION == 1 && CASACORE_PATCH_VERSION < 2))) - // In casacore < 3.1.2 thisMSAErrorHandler is a raw pointer, - // From casacore 3.1.2. it's a CountedPtr - BaselineSelectErrorHandler errorHandler(os); - MSAntennaParse::thisMSAErrorHandler = &errorHandler; -#else - // After casacore 3.5.0, the type of MSAntennaParse::thisMSAErrorHandler - // changed again from a CountedPtr to a unique_ptr, so derive the appropriate - // type: - using CasacorePointerType = decltype(MSAntennaParse::thisMSAErrorHandler); - CasacorePointerType errorHandler(new common::BaselineSelectErrorHandler(os)); - MSAntennaParse::thisMSAErrorHandler = std::move(errorHandler); -#endif - try { - // Create a table expression representing the selection. - TableExprNode node = msAntennaGramParseCommand( - anttab, a1Node, a2Node, baselineSelection, selectedAnts1, selectedAnts2, - selectedBaselines); - // Get the antenna numbers. - Table seltab = a1Node.table()(node); - casacore::Vector<int> a1 = - ScalarColumn<int>(seltab, "ANTENNA1").getColumn(); - casacore::Vector<int> a2 = - ScalarColumn<int>(seltab, "ANTENNA2").getColumn(); - int nant = anttab.nrow(); - Matrix<bool> bl(nant, nant, false); - for (unsigned int i = 0; i < a1.size(); ++i) { - bl(a1[i], a2[i]) = true; - bl(a2[i], a1[i]) = true; - } - MSAntennaParse::thisMSAErrorHandler = curHandler; - return bl; - } catch (const std::exception&) { - MSAntennaParse::thisMSAErrorHandler = curHandler; - throw; - } -} - -BaselineSelectErrorHandler::~BaselineSelectErrorHandler() {} - -void BaselineSelectErrorHandler::reportError(const char* token, - const casacore::String message) { - itsStream << message << token << '\n'; -} - -} // namespace common -} // namespace dp3 diff --git a/common/BaselineSelect.h b/common/BaselineSelect.h deleted file mode 100644 index 01de8306a..000000000 --- a/common/BaselineSelect.h +++ /dev/null @@ -1,77 +0,0 @@ -// BaselineSelect.h: Convert MSSelection baseline string to a Matrix -// -// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy) -// SPDX-License-Identifier: GPL-3.0-or-later - -/// @file -/// Convert MSSelection baseline string to a Matrix -/// @author Ger van Diepen (diepen AT astron nl) - -#ifndef MS_BASELINESELECT_H -#define MS_BASELINESELECT_H - -#include <casacore/ms/MSSel/MSSelectionErrorHandler.h> -#include <casacore/casa/Arrays/Array.h> - -#include <ostream> - -namespace casacore { -class Table; -class TableExprNode; -class MPosition; -} // namespace casacore - -namespace dp3 { -namespace common { - -/// @ingroup MS -/// @brief Convert MSSelection baseline string to a Matrix -/// @{ - -/// Class with a static function to convert a casacore MSSelection baseline -/// string to a Matrix<Bool> telling which baselines are selected. - -class BaselineSelect { - public: - /// Parse an MSSelection baseline string and create a Matrix telling - /// which baselines are selected. - /// @param names Name for each station/antenna. - /// @param positions Position for each station/antenna. - /// @param antenna1 For each baseline, the index of the first antenna. - /// @param antenna2 For each baseline, the index of the second antenna. - /// @param baseline_selection Selection string in MSSelection format. - /// @param os Possible messages from the parser are written to this stream. - /// @return An n_stations x n_stations Matrix, which holds true for selected - // baselines and false for the other baselines. - static casacore::Matrix<bool> convert( - const std::vector<std::string>& names, - const std::vector<casacore::MPosition>& positions, - const std::vector<int>& antenna1, const std::vector<int>& antenna2, - const std::string& baseline_selection, std::ostream&); - - private: - static casacore::Matrix<bool> convert(casacore::Table& anttab, - casacore::TableExprNode& a1, - casacore::TableExprNode& a2, - const string& baselineSelection, - std::ostream& os); -}; - -/// @brief This class handles an error from the Casacore's MSAntennaParse. -/// It adds the message to the message list of the parent BaselineSelect. -class BaselineSelectErrorHandler : public casacore::MSSelectionErrorHandler { - public: - BaselineSelectErrorHandler(std::ostream& os) : itsStream(os) {} - ~BaselineSelectErrorHandler() override; - void reportError(const char* token, const casacore::String message) override; - - private: - std::ostream& itsStream; -}; - -/// @} - -} // namespace common -} // namespace dp3 - -#endif diff --git a/steps/MSBDAReader.cc b/steps/MSBDAReader.cc index bcca20517..936d7c4bc 100644 --- a/steps/MSBDAReader.cc +++ b/steps/MSBDAReader.cc @@ -12,7 +12,6 @@ #include "../base/MS.h" #include "../common/ParameterSet.h" -#include "../common/BaselineSelect.h" #include <casacore/casa/OS/HostInfo.h> #include <casacore/tables/Tables/TableRecord.h> diff --git a/steps/MSReader.cc b/steps/MSReader.cc index e4e2775b2..f4a4f7c91 100644 --- a/steps/MSReader.cc +++ b/steps/MSReader.cc @@ -35,12 +35,11 @@ #include <casacore/casa/Quanta/MVTime.h> #include <casacore/casa/OS/Conversion.h> -#include "../base/MS.h" +#include <aocommon/logger.h> +#include "../base/BaselineSelection.h" +#include "../base/MS.h" #include "../common/ParameterSet.h" -#include "../common/BaselineSelect.h" - -#include <aocommon/logger.h> using casacore::ArrayColumn; using casacore::ArrayMeasColumn; @@ -132,45 +131,21 @@ MSReader::MSReader(const casacore::MeasurementSet& ms, MSSelection select; // Overwrite the error handler to ignore errors for unknown antennas. - // Borrowed from BaselineSelect.cc - std::ostringstream os; - auto curHandler = MSAntennaParse::thisMSAErrorHandler; -#if CASACORE_MAJOR_VERSION < 3 || \ - (CASACORE_MAJOR_VERSION == 3 && \ - (CASACORE_MINOR_VERSION == 0 || \ - (CASACORE_MINOR_VERSION == 1 && CASACORE_PATCH_VERSION < 2))) - // In casacore < 3.1.2 thisMSAErrorHandler is a raw pointer, - // From casacore 3.1.2. it's a CountedPtr - common::BaselineSelectErrorHandler errorHandler(os); - MSAntennaParse::thisMSAErrorHandler = &errorHandler; -#else - // After casacore 3.5.0, the type of MSAntennaParse::thisMSAErrorHandler - // changed again from a CountedPtr to a unique_ptr, so derive the - // appropriate type: - using CasacorePointerType = decltype(MSAntennaParse::thisMSAErrorHandler); - CasacorePointerType errorHandler( - new common::BaselineSelectErrorHandler(os)); - MSAntennaParse::thisMSAErrorHandler = std::move(errorHandler); -#endif + // First construct MSSelection, because it resets the error handler. + dp3::base::LogAntennaParseErrors ignore_unknown_antennas; // Set given selection strings. - try { - select.setAntennaExpr(itsSelBL); - // Create a table expression for an MS representing the selection. - MeasurementSet ms(itsSelMS); - TableExprNode node = select.toTableExprNode(&ms); - Table subset = itsSelMS(node); - // If not all is selected, use the selection. - if (subset.nrow() < itsSelMS.nrow()) { - if (subset.nrow() <= 0) - throw std::runtime_error("Baselines " + itsSelBL + "not found in " + - msName()); - itsSelMS = subset; - } - MSAntennaParse::thisMSAErrorHandler = curHandler; - } catch (const std::exception&) { - MSAntennaParse::thisMSAErrorHandler = curHandler; - throw; + select.setAntennaExpr(itsSelBL); + // Create a table expression for an MS representing the selection. + MeasurementSet ms(itsSelMS); + TableExprNode node = select.toTableExprNode(&ms); + Table subset = itsSelMS(node); + // If not all is selected, use the selection. + if (subset.nrow() < itsSelMS.nrow()) { + if (subset.nrow() <= 0) + throw std::runtime_error("Baselines " + itsSelBL + "not found in " + + msName()); + itsSelMS = subset; } } // Prepare the MS access and get time info. -- GitLab