From 8b026a0796c8c2b2938dd71cf07a73b97a3af072 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <veenboer@astron.nl> Date: Mon, 13 Jun 2022 13:19:22 +0000 Subject: [PATCH 01/61] PADRE-1: Add AntennaFlagger step The newly added AntennaFlagger step implements the preflagger.baselines functionality of the PreFlagger step, with a better performing alternative implementation based on regular expressions. tAntennaFlagger provides a basic test for this new step. --- CMakeLists.txt | 2 + base/DP3.cc | 4 + common/baseline_utils.h | 79 ++++++++++++ steps/AntennaFlagger.cc | 189 +++++++++++++++++++++++++++++ steps/AntennaFlagger.h | 40 ++++++ steps/test/unit/tAntennaFlagger.cc | 182 +++++++++++++++++++++++++++ 6 files changed, 496 insertions(+) create mode 100644 common/baseline_utils.h create mode 100644 steps/AntennaFlagger.cc create mode 100644 steps/AntennaFlagger.h create mode 100644 steps/test/unit/tAntennaFlagger.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index c7562910b..6d1dc9abb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -384,6 +384,7 @@ add_library( base/SubtractMixed.cc base/SubtractNew.cc base/UVWCalculator.cc + steps/AntennaFlagger.cc steps/AOFlaggerStep.cc steps/ApplyCal.cc steps/Averager.cc @@ -579,6 +580,7 @@ if(BUILD_TESTING) ddecal/test/unit/tSolvers.cc steps/test/unit/mock/MockInput.cc steps/test/unit/mock/MockStep.cc + steps/test/unit/tAntennaFlagger.cc steps/test/unit/tAOFlaggerStep.cc steps/test/unit/tApplyCal.cc steps/test/unit/tApplyCalH5.cc diff --git a/base/DP3.cc b/base/DP3.cc index 7511c8372..d48ec2ad9 100644 --- a/base/DP3.cc +++ b/base/DP3.cc @@ -13,6 +13,7 @@ #include <boost/algorithm/string.hpp> +#include "../steps/AntennaFlagger.h" #include "../steps/AOFlaggerStep.h" #include "../steps/ApplyBeam.h" #include "../steps/ApplyCal.h" @@ -52,6 +53,7 @@ #include "../common/Timer.h" #include "../common/StreamUtil.h" +#include "../steps/AntennaFlagger.h" #include <casacore/casa/OS/Path.h> #include <casacore/casa/OS/DirectoryIterator.h> @@ -161,6 +163,8 @@ static std::shared_ptr<Step> makeSingleStep(const std::string& type, step = std::make_shared<steps::MedFlagger>(inputStep, parset, prefix); } else if (type == "preflagger" || type == "preflag") { step = std::make_shared<steps::PreFlagger>(inputStep, parset, prefix); + } else if (type == "antennaflagger") { + step = std::make_shared<steps::AntennaFlagger>(inputStep, parset, prefix); } else if (type == "uvwflagger" || type == "uvwflag") { step = std::make_shared<steps::UVWFlagger>(inputStep, parset, prefix, inputType); diff --git a/common/baseline_utils.h b/common/baseline_utils.h new file mode 100644 index 000000000..70f320b4b --- /dev/null +++ b/common/baseline_utils.h @@ -0,0 +1,79 @@ +/* + * File taken from Aartfaac2ms repository: + * https://git.astron.nl/RD/aartfaac-tools/-/blob/master/common/baseline_utils.h + */ + +#ifndef LOFAR_COMMON_BASELINE_UTILS_H_ +#define LOFAR_COMMON_BASELINE_UTILS_H_ + +#include <cstdint> +#include <utility> + +namespace dp3 { +namespace common { + +/* + * Computes the number of baselines, including auto-correlations. + */ +inline int ComputeNBaselines(int nAntennas) { + return nAntennas * (nAntennas - 1) / 2 + nAntennas; +} + +/* + * Computes the baseline index corresponding to antenna a and b, given the + * number of antennas. It assumes that auto-correlations of antennas are + * included. + */ +inline int ComputeBaselineIndex(int a, int b, int nAntennas) { + return (a * nAntennas) + b - a * (a + 1) / 2; +} + +/* + * Computes the antenna pair (baseline) given the baseline index. It assumes + * that auto-correlations of antennas are included. + */ +inline std::pair<int, int> ComputeBaseline(int baselineIndex, int nAntennas) { + const int nBaselines = ComputeNBaselines(nAntennas); + // Compute first antenna for baselineIndex + const int antenna1 = + nAntennas - (-1 + std::sqrt(8 * (nBaselines - baselineIndex))) / 2; + + // Compute the number of baselines (antenna1,*) + const int nBaselinesAntenna1 = + ((nAntennas * 2 + 1) * antenna1 - (antenna1 * antenna1)) / 2; + + // Compute the second antenna for baselineIndex + const int antenna2 = baselineIndex - nBaselinesAntenna1 + antenna1; + + return {antenna1, antenna2}; +} + +/* + * Computes the list of baselines for a given antenna. It assumes that + * auto-correlations are included. + */ +inline std::vector<int> ComputeBaselineList(int antenna, int nAntennas) { + std::vector<int> indices(nAntennas); + + for (int i = antenna; i < nAntennas; i++) { + // Compute the number of baselines up to antenna i + int m = int((i * (i + 1)) / 2); + + if (i == antenna) { + // Add all baselines (j, antenna) + for (int j = 0; j < antenna + 1; j++) { + indices[j] = m + j; + } + } else { + // Add baseline (antenna, i) + indices[i] = m + antenna; + } + } + + return indices; +} + +} // namespace common +} // namespace dp3 + +#endif // LOFAR_COMMON_BASELINE_UTILS_H_ diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc new file mode 100644 index 000000000..0011914be --- /dev/null +++ b/steps/AntennaFlagger.cc @@ -0,0 +1,189 @@ + +#include "AntennaFlagger.h" +#include "../common/ParameterSet.h" +#include "../common/baseline_utils.h" + +#include <boost/algorithm/string.hpp> +#include <regex> + +namespace { +/** + * Converts an input string with the format defined in preflagger's baseline + * selection formats. + * This can is further explained in + * https://www.astron.nl/citt/DP3/steps/Description%20of%20baseline%20selection%20parameters.html + * + * Currently, the following selections are implemented; + * - Single antennas (e.g. 1,3,15,4) + * - Antenna ranges (e.g. 1~15) + * - Autocorrelations for a selection (e.g. 1~15&&& does all + * autocorrelations for the antennas in a selection) + * Multiple selections can be separated with a semicolon (;) + * @param selection the input string to-be parsed. + * @param n_antennas the total number of antennas + * @return a vector of pairs in the format (ant1, ant2) + */ +std::vector<std::pair<int, int>> convertSelection( + const std::string& selectionString, int nAntennas) { + // The list of all baselines to be flagged + std::vector<std::pair<int, int>> selection; + + // Split selection string on antenna-sets + std::vector<std::string> antenna_sets; + boost::split(antenna_sets, selectionString, boost::is_any_of(";")); + + for (std::string& antenna_set : antenna_sets) { + // Split on antennas within the antenna-set. + std::vector<std::string> antenna_strings; + boost::split(antenna_strings, antenna_set, boost::is_any_of(",")); + + // Check if correlation mode for this set of antennas. + bool auto_correlation = false; + bool cross_correlation = false; + std::smatch matches; + if (std::regex_search(antenna_set, matches, std::regex("&&&"))) { + auto_correlation = true; + } else if (std::regex_search(antenna_set, matches, std::regex("&&"))) { + auto_correlation = true; + cross_correlation = true; + } else if (std::regex_search(antenna_set, matches, std::regex("&"))) { + cross_correlation = true; + } else { // Nothing specified, same as && + auto_correlation = true; + cross_correlation = true; + } + + // Convert all antenna numbers from input string to a vector of ints. + std::vector<int> antennas; + for (const std::string& antenna : antenna_strings) { + // This regex searches for a range, low~high (e.g. 0~20). + std::regex range("(\\d+)~(\\d+)"); + std::smatch matches; + + try { + if (std::regex_search(antenna, matches, range)) { + // A range like 0~20 is in the string. + int low = std::stoi(matches[1]); + int high = std::stoi(matches[2]) + 1; // High is inclusive + for (; low < high; low++) { + antennas.push_back(low); + } + } else { + antennas.push_back(std::stoi(antenna)); + } + } catch (std::invalid_argument& e) { + std::cerr << "Could not parse: " << antenna << std::endl; + continue; + } + } + + // Iterate antennas and create tuple baselines and insert into vector. + if (auto_correlation) { + // Ignore autocorrelated antenna indices. + for (int antenna : antennas) { + selection.emplace_back(antenna, antenna); + } + } + + // Ignore both baseline combinations for the antennas in the vector. + if (cross_correlation) { + for (int antenna1 : antennas) { + for (int antenna2 = 0; antenna2 < antenna1; antenna2++) { + selection.emplace_back(antenna1, antenna2); + } + } + } + } + + return selection; +} + +casacore::Matrix<bool> convert(const std::string& selectionString, + uint nAntennas) { + assert(nAntennas != 0); + casacore::Matrix<bool> bl(nAntennas, nAntennas, false); + + std::vector<std::pair<int, int>> selection = + convertSelection(selectionString, int(nAntennas)); + + for (std::pair<int, int> p : selection) { + bl(p.first, p.second) = true; + bl(p.second, p.first) = true; + } + return bl; +} +} // namespace + +namespace dp3 { +namespace steps { + +void AntennaFlagger::finish() { getNextStep()->finish(); } + +void AntennaFlagger::show(std::ostream& ostream) const { + ostream << "AntennaFlagger " << name_; + ostream << "\n selection: " << selectionString_; +} + +AntennaFlagger::AntennaFlagger(InputStep* input, + const common::ParameterSet& parset, + const std::string& prefix) + : name_(prefix) { + initializationTimer_.start(); + + // Parse arguments + std::string selectionString = + parset.getString(prefix + "selection", string()); + + if (selectionString.empty()) { + throw std::runtime_error("Selection needs to be set for AntennaFlagger."); + } + + // Parse our selection string into a matrix of booleans. + n_antennas_ = input->getInfo().nantenna(); + selectionString_ = selectionString; + selection_ = convert(selectionString_, n_antennas_); + + initializationTimer_.stop(); +} + +bool AntennaFlagger::process(const base::DPBuffer& inBuffer) { + computationTimer_.start(); + + buffer_.copy(inBuffer); + + // Flags are in format: corr,chan,baseline. + casacore::Cube<bool>& flags = buffer_.getFlags(); + + for (size_t correlation = 0; correlation < flags.nrow(); correlation++) { + for (size_t channel = 0; channel < flags.ncolumn(); channel++) { + for (size_t baseline = 0; baseline < flags.nplane(); baseline++) { + const std::pair<size_t, size_t> antennas = + common::ComputeBaseline(baseline, n_antennas_); + flags(correlation, channel, baseline) = + selection_(antennas.first, antennas.second); + } + } + } + + computationTimer_.stop(); + + getNextStep()->process(buffer_); + + return true; +} + +void AntennaFlagger::showTimings(std::ostream& ostream, double duration) const { + double initTime = initializationTimer_.getElapsed(); + double computeTime = computationTimer_.getElapsed(); + double totTime = initTime + computeTime; + + ostream << " "; + base::FlagCounter::showPerc1(ostream, totTime, duration); + ostream << " AntennaFlagger " << name_ << "\n "; + base::FlagCounter::showPerc1(ostream, initTime, totTime); + ostream << " of it spent in initialization.\n "; + base::FlagCounter::showPerc1(ostream, computeTime, totTime); + ostream << " of it spent in computation.\n"; +} +} // namespace steps +} // namespace dp3 \ No newline at end of file diff --git a/steps/AntennaFlagger.h b/steps/AntennaFlagger.h new file mode 100644 index 000000000..b5e85ce4a --- /dev/null +++ b/steps/AntennaFlagger.h @@ -0,0 +1,40 @@ +#ifndef DPPP_ANTENNAFLAGGER_H +#define DPPP_ANTENNAFLAGGER_H + +#include "InputStep.h" + +#include "../base/DPBuffer.h" +#include "../base/BaselineSelection.h" + +#include <casacore/measures/Measures/MDirection.h> + +namespace dp3 { +namespace common { +class ParameterSet; +class ParameterValue; +} // namespace common + +namespace steps { +class AntennaFlagger final : public Step { + public: + AntennaFlagger(InputStep* input, const common::ParameterSet& parset, + const string& prefix); + void finish() override; + void show(std::ostream& ostream) const override; + bool process(const base::DPBuffer& buffer) override; + void showTimings(std::ostream& ostream, double duration) const override; + + private: + string name_; + base::DPBuffer buffer_; + size_t n_antennas_; + std::string selectionString_; + casacore::Matrix<bool> selection_; + common::NSTimer initializationTimer_; + common::NSTimer computationTimer_; +}; + +} // namespace steps +} // namespace dp3 + +#endif diff --git a/steps/test/unit/tAntennaFlagger.cc b/steps/test/unit/tAntennaFlagger.cc new file mode 100644 index 000000000..76caa0d6d --- /dev/null +++ b/steps/test/unit/tAntennaFlagger.cc @@ -0,0 +1,182 @@ +// tAntennaFlagger.cc: Test program for class AntennaFlagger +// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy) +// SPDX-License-Identifier: GPL-3.0-or-later + +#include <casacore/casa/Arrays/ArrayMath.h> +#include <casacore/casa/Arrays/ArrayLogical.h> + +#include <boost/test/unit_test.hpp> + +#include "tStepCommon.h" +#include "../../AntennaFlagger.h" +#include "../../Counter.h" +#include "../../InputStep.h" +#include "../../../base/DPBuffer.h" +#include "../../../base/DPInfo.h" +#include "../../../common/ParameterSet.h" +#include "../../../common/StringTools.h" + +using dp3::base::DPBuffer; +using dp3::base::DPInfo; +using dp3::common::ParameterSet; +using dp3::steps::AntennaFlagger; +using dp3::steps::Counter; +using dp3::steps::InputStep; +using dp3::steps::Step; + +BOOST_AUTO_TEST_SUITE(antennaflagger) + +// Simple class to generate input arrays. +// It can only set all flags to true or all to false. +// Weights are always 1. +// It can be used with different nr of times, channels, etc. +class TestInput final : public InputStep { + public: + TestInput(int ntime, int nbl, int nchan, int ncorr, bool flag) + : count_(0), + n_times_(ntime), + n_bl_(nbl), + n_chan_(nchan), + n_corr__(ncorr), + flag_(flag) { + // Define start time 0.5 (= 3 - 0.5*5) and time interval 5. + info().init(ncorr, 0, nchan, ntime, 0.5, 5., string(), string()); + + // Fill the baselines; use 12 stations with 48 antennas each. + constexpr int nreceivers = 12 * 48; + std::vector<int> ant1(nbl); + std::vector<int> ant2(nbl); + int st1 = 0; + int st2 = 0; + for (int i = 0; i < nbl; ++i) { + ant1[i] = st1; + ant2[i] = st2; + if (++st2 == nreceivers) { + st2 = 0; + if (++st1 == nreceivers) { + st1 = 0; + } + } + } + + // Atenna names and positions are not used by the AntennaFlagger, + // but the arrays need to have the correct dimensions. + std::vector<std::string> antNames(nreceivers); + std::vector<casacore::MPosition> antPos(nreceivers); + std::vector<double> antDiam(nreceivers, 70.); + info().set(antNames, antDiam, antPos, ant1, ant2); + } + + private: + bool process(const DPBuffer&) override { + // Stop when all times are done. + if (count_ == n_times_) { + return false; + } + casacore::Cube<casacore::Complex> data(n_corr__, n_chan_, n_bl_); + for (int i = 0; i < int(data.size()); ++i) { + data.data()[i] = casacore::Complex(i + count_ * 10, i - 10 + count_ * 6); + } + casacore::Matrix<double> uvw(3, n_bl_); + for (int i = 0; i < n_bl_; ++i) { + uvw(0, i) = 1 + count_ + i; + uvw(1, i) = 2 + count_ + i; + uvw(2, i) = 3 + count_ + i; + } + DPBuffer buf; + buf.setTime(count_ * 5 + 3); // same interval as in updateAveragInfo + buf.setData(data); + buf.setUVW(uvw); + casacore::Cube<float> weights(data.shape()); + weights = 1.0; + buf.setWeights(weights); + casacore::Cube<bool> flags(data.shape()); + flags = flag_; + buf.setFlags(flags); + // The fullRes flags are a copy of the XX flags, but differently shaped. + // They are not averaged, thus only 1 time per row. + casacore::Cube<bool> fullResFlags(n_chan_, 1, n_bl_); + fullResFlags = flag_; + buf.setFullResFlags(fullResFlags); + getNextStep()->process(buf); + ++count_; + return true; + } + + void finish() override { getNextStep()->finish(); } + void show(std::ostream&) const override{}; + void updateInfo(const DPInfo&) override {} + + int count_; + int n_times_; + int n_bl_; + int n_chan_; + int n_corr__; + bool flag_; +}; + +// Class to check result of flagged, unaveraged TestInput run by test1. +class TestOutput : public Step { + public: + TestOutput(int ntime, int nbl, int nchan, int ncorr) + : count_(0), + n_times_(ntime), + n_bl_(nbl), + n_chan_(nchan), + n_corr_(ncorr) {} + + private: + virtual bool process(const DPBuffer& buf) { + casacore::Cube<bool> result(n_corr_, n_chan_, n_bl_); + for (int i = 0; i < n_bl_; ++i) { + casacore::Cube<bool> test = buf.getFlags(); + + if (i == 4 || i == 10 || i == 11) { + for (int j = 0; j < n_chan_; ++j) { + for (int k = 0; k < n_corr_; ++k) { + result(k, j, i) = true; + } + } + } + } + BOOST_CHECK(allEQ(buf.getFlags(), result)); + count_++; + return true; + } + + virtual void finish() {} + virtual void show(std::ostream&) const {} + virtual void updateInfo(const DPInfo& infoIn) { + info() = infoIn; + BOOST_CHECK_EQUAL(int(infoIn.origNChan()), n_chan_); + BOOST_CHECK_EQUAL(int(infoIn.nchan()), n_chan_); + BOOST_CHECK_EQUAL(int(infoIn.ntime()), n_times_); + BOOST_CHECK_EQUAL(infoIn.timeInterval(), 5); + BOOST_CHECK_EQUAL(int(infoIn.nchanAvg()), 1); + BOOST_CHECK_EQUAL(int(infoIn.ntimeAvg()), 1); + } + + int count_; + int n_times_; + int n_bl_; + int n_chan_; + int n_corr_; + bool flag_; + bool clear_; +}; + +// Test flagging a few antennae. +void test1(int ntime, int nbl, int nchan, int ncorr, bool flag) { + // Create the steps. + TestInput* in = new TestInput(ntime, nbl, nchan, ncorr, flag); + Step::ShPtr step1(in); + ParameterSet parset; + parset.add("selection", "4,10,11"); + Step::ShPtr step2 = std::make_shared<AntennaFlagger>(in, parset, ""); + Step::ShPtr step3 = std::make_shared<TestOutput>(ntime, nbl, nchan, ncorr); + dp3::steps::test::Execute({step1, step2, step3}); +} + +BOOST_AUTO_TEST_CASE(test_1) { test1(10, 576, 8, 4, false); } + +BOOST_AUTO_TEST_SUITE_END() -- GitLab From d5bd680f3863cfd139f79c83758e82e6b0fd39d5 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <veenboer@astron.nl> Date: Wed, 6 Jul 2022 12:27:36 +0200 Subject: [PATCH 02/61] PADRE-2: Add alternative AntennaFlagger We integrated the flagging functionality of the A12 pipeline, which consists of separate calls (aoquality, aoflagger, a Python script and DP3) into the AntennaFlagger step. This approach is not only much simpler, but also much faster. Antennaflagger is extended with a 'flagging' parameter to flag the input data. --- .gitmodules | 3 +- CMakeLists.txt | 14 +- antennaflagger/AntennaFlagger.cc | 411 +++++++++++++++++++++++++++++ antennaflagger/AntennaFlagger.h | 77 ++++++ antennaflagger/CMakeLists.txt | 1 + antennaflagger/statistics.h | 228 ++++++++++++++++ common/baseline_utils.h | 125 +++++++-- common/test/unit/tBaselineUtils.cc | 118 +++++++++ docs/schemas/AntennaFlagger.yml | 24 ++ steps/AntennaFlagger.cc | 133 ++++++++-- steps/AntennaFlagger.h | 42 +-- steps/test/unit/tAntennaFlagger.cc | 12 +- 12 files changed, 1110 insertions(+), 78 deletions(-) create mode 100644 antennaflagger/AntennaFlagger.cc create mode 100644 antennaflagger/AntennaFlagger.h create mode 100644 antennaflagger/CMakeLists.txt create mode 100644 antennaflagger/statistics.h create mode 100644 common/test/unit/tBaselineUtils.cc create mode 100644 docs/schemas/AntennaFlagger.yml diff --git a/.gitmodules b/.gitmodules index 47a248d6e..373086ebe 100644 --- a/.gitmodules +++ b/.gitmodules @@ -7,7 +7,6 @@ [submodule "external/pybind11"] path = external/pybind11 url = https://github.com/pybind/pybind11.git - [submodule "external/schaap-packaging"] path = external/schaap-packaging - url = https://git.astron.nl/RD/schaap-packaging.git + url = https://git.astron.nl/RD/schaap-packaging.git \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 6d1dc9abb..a2d539708 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -322,6 +322,9 @@ add_library( common/VdsMaker.cc common/VdsPartDesc.cc) +add_library(AntennaFlagger OBJECT antennaflagger/AntennaFlagger.cc) +set_property(TARGET AntennaFlagger PROPERTY POSITION_INDEPENDENT_CODE ON) + add_library( DDECal OBJECT steps/BdaDdeCal.cc @@ -458,9 +461,13 @@ set_property(TARGET Blob Common DDECal DP3_OBJ ParmDB PythonDP3 PROPERTY POSITION_INDEPENDENT_CODE ON) set(DP3_OBJECTS - $<TARGET_OBJECTS:Blob> $<TARGET_OBJECTS:Common> $<TARGET_OBJECTS:DDECal> - $<TARGET_OBJECTS:DP3_OBJ> $<TARGET_OBJECTS:ParmDB> - $<TARGET_OBJECTS:PythonDP3>) + $<TARGET_OBJECTS:Blob> + $<TARGET_OBJECTS:Common> + $<TARGET_OBJECTS:DDECal> + $<TARGET_OBJECTS:DP3_OBJ> + $<TARGET_OBJECTS:ParmDB> + $<TARGET_OBJECTS:PythonDP3> + $<TARGET_OBJECTS:AntennaFlagger>) set(LIBDIRAC_PREFIX "" @@ -617,6 +624,7 @@ if(BUILD_TESTING) list( APPEND TEST_FILENAMES + common/test/unit/tBaselineUtils.cc common/test/unit/tMemory.cc common/test/unit/tStringTools.cc base/test/unit/tDPInfo.cc diff --git a/antennaflagger/AntennaFlagger.cc b/antennaflagger/AntennaFlagger.cc new file mode 100644 index 000000000..54fd5809d --- /dev/null +++ b/antennaflagger/AntennaFlagger.cc @@ -0,0 +1,411 @@ +// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy) +// SPDX-License-Identifier: GPL-3.0-or-later + +#if defined(DEBUG) +#include <cassert> +#endif +#include <cstddef> +#include <iostream> +#include <set> + +#include "AntennaFlagger.h" +#include "./statistics.h" + +// Non-member helper functions +namespace { +inline size_t data_index(int n_channels, int n_correlations, int bl, int chan, + int cor) { + // data[n_baselines][n_channels][n_correlations] + return (bl * n_channels + chan) * n_correlations + cor; +} + +/* + * Helper functions for AaartfaacFlagger::find_bad_antennas + */ +template <typename InType, typename OutType> +void compute_stats_per_baseline( + size_t n_baselines, size_t n_channels, size_t n_correlations, + const std::vector<std::complex<InType>>& data, + std::vector<std::complex<OutType>>* stats_std, + std::vector<std::complex<OutType>>* stats_sump2) { + // - in: data[n_baselines][n_channels][n_correlations] + // - out: stats[n_correlations][n_baselines] + // Note that the baseline and correlation indices are swapped between + // input and output. This is also different from the Python reference, + // and is to avoid strided access later on where the stats are used. + +#if defined(DEBUG) + size_t n_in = n_baselines * n_channels * n_correlations; + assert(data.size() == n_in); +#endif + + std::vector<OutType> result(n_baselines * n_correlations); + size_t n_out = n_correlations * n_baselines; + stats_std->resize(n_out); + stats_sump2->resize(n_out); + + for (size_t bl = 0; bl < n_baselines; bl++) { + for (size_t cor = 0; cor < n_correlations; cor++) { + // Compute sum, sump2 + OutType sum_real = 0; + OutType sum_imag = 0; + OutType sump2_real = 0; + OutType sump2_imag = 0; + for (size_t chan = 0; chan < n_channels; chan++) { + size_t idx = data_index(n_channels, n_correlations, bl, chan, cor); + std::complex<InType> value = data[idx]; + sum_real += value.real(); + sum_imag += value.imag(); + sump2_real += value.real() * value.real(); + sump2_imag += value.imag() * value.imag(); + } + + // Compute mean + OutType mean_real = sum_real / n_channels; + OutType mean_imag = sum_imag / n_channels; + + // Compute standard deviation + OutType sum2_real = 0; + OutType sum2_imag = 0; + for (size_t chan = 0; chan < n_channels; chan++) { + size_t idx = data_index(n_channels, n_correlations, bl, chan, cor); + std::complex<InType> value = data[idx]; + OutType diff_real = mean_real - value.real(); + OutType diff_imag = mean_imag - value.imag(); + sum2_real += diff_real * diff_real; + sum2_imag += diff_imag * diff_imag; + } + OutType std_real = std::sqrt(sum2_real / n_channels); + OutType std_imag = std::sqrt(sum2_imag / n_channels); + + // Store result + size_t idx = cor * n_baselines + bl; + (*stats_std)[idx] = std::complex<OutType>(std_real, std_imag); + (*stats_sump2)[idx] = std::complex<OutType>(sump2_real, sump2_imag); + } + } +} // end compute_stats_per_baseline + +template <typename InType, typename OutType> +void compute_stats_per_antenna( + int n_receivers, int n_baselines, int n_channels, int n_correlations, + const std::vector<std::complex<InType>>& data, + std::vector<std::complex<OutType>>* stats_std_grouped, + std::vector<std::complex<OutType>>* stats_sump2_grouped) { + // Compute stats for every baseline and correlation: + // - in: data[n_baselines][n_channels][n_correlations] + // - out: stats_std[n_correlations][n_baselines] + // - out: stats_sump2[n_correlations][n_baselines] +#if defined(DEBUG) + assert(data.size() == size_t(n_baselines * n_channels * n_correlations)); +#endif + std::vector<std::complex<OutType>> stats_std; + std::vector<std::complex<OutType>> stats_sump2; + compute_stats_per_baseline(n_baselines, n_channels, n_correlations, data, + &stats_std, &stats_sump2); +#if defined(DEBUG) + assert(stats_std.size() == size_t(n_correlations * n_baselines)); + assert(stats_sump2.size() == size_t(n_correlations * n_baselines)); +#endif + + // Group stats by antenna: + // - in: stats_std[n_correlations][n_baselines] + // - in: stats_sump2[n_correlations][n_baselines] + // - out: stats_std_grouped[n_correlations][n_receivers] + // - out: stats_sump2_grouped[n_correlations][n_receivers] + stats_std_grouped->resize(n_receivers * n_correlations); + stats_sump2_grouped->resize(n_receivers * n_correlations); + + for (int antenna = 0; antenna < n_receivers; antenna++) { + // Find all baselines with this antenna + std::vector<int> indices = dp3::common::ComputeBaselineList( + antenna, n_receivers, dp3::common::BaselineOrder::ROW_MAJOR); + + // Sum statistics over these baselines + for (int cor = 0; cor < n_correlations; cor++) { + for (const int& bl : indices) { + size_t idx_src = cor * n_baselines + bl; + size_t idx_dst = cor * n_receivers + antenna; + (*stats_std_grouped)[idx_dst] += stats_std[idx_src]; + (*stats_sump2_grouped)[idx_dst] += stats_sump2[idx_src]; + } + } + } +} // end compute_stats_per_antenna + +/* + * Helper functions for AaartfaacFlagger::find_bad_stations + */ +template <typename T> +void find_bad_antennas_(int n_stations, int n_receivers_per_station, + int n_baselines, int n_correlations, float sigma, + int maxiters, const std::vector<std::complex<T>>& stats, + std::vector<int>* all_bad_stations) { + int n_receivers = n_stations * n_receivers_per_station; + // - in: stats[n_correlations][n_receivers] + // - out: all_bad_stations[?] +#if defined(DEBUG) + assert(n_correlations == static_cast<int>(4)); + assert(stats.size() == static_cast<size_t>(n_correlations * n_receivers)); +#endif + + // Find bad stations by looking at the statistics (per correlation) for all + // antennas that belong to one station. The results are stored in a list of + // lists such that stations can be processed independently. + std::vector<std::vector<int>> bad_stations(n_stations); + + // Check XY, YX and YY + for (int cor = 1; cor < n_correlations; cor++) { + // The check sum is the list of statistics for the current correlation + // - in: T2 check_sum[n_correlations][n_receivers] + const std::complex<T>* check_sum = &stats[cor * n_receivers]; + + for (int station = 0; station < n_stations; station++) { + // Make a list of all antennas for the current station, + // and select the corresponding values from check_sum. + std::vector<int> station_indices(n_receivers_per_station); + std::vector<T> check_sum_station(n_receivers_per_station * 2); + for (int i = 0; i < n_receivers_per_station; i++) { + int station_index = station * n_receivers_per_station + i; + station_indices[i] = station_index; + T real = std::log(check_sum[station_index].real()); + T imag = std::log(check_sum[station_index].imag()); + check_sum_station[0 * n_receivers_per_station + i] = real; + check_sum_station[1 * n_receivers_per_station + i] = imag; + } + + // Find outliers + int n_flags = n_receivers_per_station * 2; + std::vector<bool> flags = + find_outliers(sigma, maxiters, check_sum_station); +#if defined(DEBUG) + assert(flags.size() == size_t(n_flags)); +#endif + + // Count number of outliers + int n_outliers = 0; + for (int i = 0; i < n_flags; i++) { + n_outliers += flags[i]; + } +#if defined(DEBUG) + if (n_outliers > 0) { + std::cout << "found " << n_outliers + << " outlier antenna(s) for correlation = " << cor + << ", station = " << station << std::endl; + } +#endif + + // Add the outliers to the list of bad stations + const int n = bad_stations[station].size(); + bad_stations[station].resize(n + n_outliers); + int j = 0; + for (int i = 0; i < n_receivers_per_station * 2; i++) { + if (flags[i]) { + int station_index = + station * n_receivers_per_station + (i % n_receivers_per_station); + bad_stations[station][n + j++] = station_index; + } + } + } + } + + // Add the bad stations to the global list of bad stations + for (int station = 0; station < n_stations; station++) { + all_bad_stations->insert(all_bad_stations->end(), + bad_stations[station].begin(), + bad_stations[station].end()); + } +} // end find_bad_antennas_ + +template <typename T> +std::vector<std::complex<T>> compute_means(int n_stations, + int n_receivers_per_station, + const std::complex<T>* check_sum) { + const int n_receivers = n_stations * n_receivers_per_station; + + // - in: T2 check_sum[n_receivers] + // - out: T2 means[n_stations] + std::vector<std::complex<T>> means(n_stations); + + // Split the real and imaginary parts of the check sum + std::vector<T> check_sum_real(n_receivers); + std::vector<T> check_sum_imag(n_receivers); + for (int antenna = 0; antenna < n_receivers; antenna++) { + check_sum_real[antenna] = check_sum[antenna].real(); + check_sum_imag[antenna] = check_sum[antenna].imag(); + } + + // Compute the median for both real and imaginary parts + T median_real = compute_median(check_sum_real); + T median_imag = compute_median(check_sum_imag); + + for (int station = 0; station < n_stations; station++) { + // Make a list of all antennas for the current station, + // and select the corresponding values from check_sum + std::vector<T> check_sum_station_real(n_receivers_per_station); + std::vector<T> check_sum_station_imag(n_receivers_per_station); + for (int i = 0; i < n_receivers_per_station; i++) { + int station_index = station * n_receivers_per_station + i; + check_sum_station_real[i] = check_sum_real[station_index]; + check_sum_station_imag[i] = check_sum_imag[station_index]; + } + + // Compute the median for both the real and imaginary parts + // of the check sum for the current station + T median_station_real = compute_median(check_sum_station_real); + T median_station_imag = compute_median(check_sum_station_imag); + + // Add the mean value for the current station + means[station] = {median_station_real / median_real, + median_station_imag / median_imag}; + } + + return means; +} // end compute_means + +}; // namespace + +/* + * Member functions + */ +namespace dp3 { +namespace antennaflagger { + +void AntennaFlagger::reportRuntime() const { + std::cout << "statistics: " << computeStatisticsWatch_.getElapsed() + << ", bad_antennas: " << findbadAntennasWatch_.getElapsed() + << ", bad_stations: " << findBadStationsWatch_.getElapsed() + << std::endl; +} + +void AntennaFlagger::computeStats( + const std::vector<std::complex<double>>& data) { + computeStatisticsWatch_.start(); + compute_stats_per_antenna(n_receivers_, n_baselines_, n_channels_, + n_correlations_, data, &stats_std_, &stats_sump2_); +#if defined(DEBUG) + assert(stats_std_.size() == size_t(n_correlations_ * n_receivers_)); + assert(stats_sump2_.size() == size_t(n_correlations_ * n_receivers_)); +#endif + computeStatisticsWatch_.stop(); +} + +void AntennaFlagger::computeStats( + const std::vector<std::complex<float>>& data) { + computeStatisticsWatch_.start(); + compute_stats_per_antenna(n_receivers_, n_baselines_, n_channels_, + n_correlations_, data, &stats_std_, &stats_sump2_); +#if defined(DEBUG) + assert(stats_std_.size() == size_t(n_correlations_ * n_receivers_)); + assert(stats_sump2_.size() == size_t(n_correlations_ * n_receivers_)); +#endif + computeStatisticsWatch_.stop(); +} + +void AntennaFlagger::assertStatsComputed() const { + if (!(stats_std_.size() == size_t(n_correlations_ * n_receivers_)) || + !(stats_sump2_.size() == size_t(n_correlations_ * n_receivers_))) { + throw std::runtime_error("Call AntennaFlagger::compute_stats() first!"); + } +} + +std::vector<int> AntennaFlagger::findBadAntennas(float sigma, int maxiters, + int outlier_flag_count) { + assertStatsComputed(); + + findbadAntennasWatch_.start(); + + std::vector<int> flagged_antennas; + find_bad_antennas_(n_stations_, n_receivers_per_station, n_baselines_, + n_correlations_, sigma, maxiters, stats_sump2_, + &flagged_antennas); + find_bad_antennas_(n_stations_, n_receivers_per_station, n_baselines_, + n_correlations_, sigma, maxiters, stats_std_, + &flagged_antennas); + + // Count the number of times an antenna is flagged + std::vector<int> flagged_count(n_receivers_); + for (int antenna : flagged_antennas) { + flagged_count[antenna]++; + } + + // Make a list of antennas that are flagged multiple times + std::vector<int> flagged_antennas2; + for (int antenna = 0; antenna < n_receivers_; antenna++) { + if (flagged_count[antenna] > outlier_flag_count) { + flagged_antennas2.push_back(antenna); + } + } + + findbadAntennasWatch_.stop(); + + return flagged_antennas2; +} + +std::vector<int> AntennaFlagger::findBadStations(float sigma, int maxiters) { + assertStatsComputed(); + + findBadStationsWatch_.start(); + + // Compute means for all the check sums seperately for + // the real and imaginary part of the check sum. + // - in: STATS_TYPE2 check_sums[n_check_sums][n_receivers] + // - out: STATS_TYPE all_means[check_sums.size() * 2][n_stations] + std::array<const std::complex<STATS_TYPE>*, 4> check_sums{ + &stats_std_[0 * n_receivers_], &stats_std_[3 * n_receivers_], + &stats_sump2_[0 * n_receivers_], &stats_sump2_[3 * n_receivers_]}; + std::vector<std::vector<STATS_TYPE>> all_means(check_sums.size() * 2); + + for (std::vector<STATS_TYPE>& means : all_means) { + means.resize(n_stations_); + } + + for (size_t i = 0; i < check_sums.size(); i++) { + std::vector<std::complex<STATS_TYPE>> means = + compute_means(n_stations_, n_receivers_per_station, check_sums[i]); + for (int station = 0; station < n_stations_; station++) { + all_means[0 * check_sums.size() + i][station] = means[station].real(); + all_means[1 * check_sums.size() + i][station] = means[station].imag(); + } + } + + // Flag all stations for which all of the means are an outlier + std::vector<bool> station_flags(n_stations_, {false}); + for (size_t i = 0; i < all_means.size(); i++) { + const std::vector<STATS_TYPE>& means = all_means[i]; + + // Find outliers + std::vector<bool> flags = find_outliers(sigma, maxiters, means); +#if defined(DEBUG) + assert(flags.size() == size_t(n_stations_)); +#endif + + for (int station = 0; station < n_stations_; station++) { + station_flags[station] = + (flags[station] && (i == 0 || station_flags[station])); + } + } + + // Make a list of all the antenna indices for the flagged stations + std::vector<int> flagged_antennas; + for (int station = 0; station < n_stations_; station++) { + if (station_flags[station]) { + // Make space for the new indices + int n = flagged_antennas.size(); + flagged_antennas.resize(n + n_receivers_per_station); + + // Add the new indices + for (int i = 0; i < n_receivers_per_station; i++) { + int station_index = station * n_receivers_per_station + i; + flagged_antennas[n + i] = station_index; + } + } + } + + findBadStationsWatch_.stop(); + + return flagged_antennas; +} + +} // namespace antennaflagger +} // namespace dp3 diff --git a/antennaflagger/AntennaFlagger.h b/antennaflagger/AntennaFlagger.h new file mode 100644 index 000000000..6860a1f1a --- /dev/null +++ b/antennaflagger/AntennaFlagger.h @@ -0,0 +1,77 @@ +// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy) +// SPDX-License-Identifier: GPL-3.0-or-later + +#ifndef DPPP_ANTENNAFLAGGER_H +#define DPPP_ANTENNAFLAGGER_H + +#include <vector> +#include <complex> + +#include "../common/baseline_utils.h" +#include "../common/Timer.h" + +namespace dp3 { +namespace antennaflagger { + +/* + * AntennaFlagger uses statistics to find bad antennas for a visibility data + * set. It is based on the following Python reference code: + * https://github.com/transientskp/A12_pipeline/blob/master/pyscripts/calc_antflags.py + * and comprises four main steps: + * 1) Compute standard deviation and sump2 for every baseline and correlation + * seperately. These statistics are integrated over frequency (channels) + * 2) Group the statistics per antenna by integrating over all baselines where + * a particular antenna is involved. These statistics are stored as class + * members. + * 3) Query these statistics to find outlier antennas + * 4) Query these statistics to find outlier stations + */ +class AntennaFlagger { + public: + AntennaFlagger(int n_stations, int n_receivers_per_station, int n_channels, + int n_correlations) + : n_stations_(n_stations), + n_receivers_per_station(n_receivers_per_station), + n_channels_(n_channels), + n_correlations_(n_correlations), + n_receivers_(n_stations * n_receivers_per_station), + n_baselines_(dp3::common::ComputeNBaselines(n_receivers_)) {} + + void computeStats(const std::vector<std::complex<double>>& data); + void computeStats(const std::vector<std::complex<float>>& data); + + std::vector<int> findBadAntennas(float sigma = 3, int maxiters = 5, + int outlier_flag_count = 2); + std::vector<int> findBadStations(float sigma = 2.5, int maxiters = 5); + + void reportRuntime() const; + + private: + // Helper function + void assertStatsComputed() const; + + // Constants + const int n_stations_; + const int n_receivers_per_station; + const int n_channels_; + const int n_correlations_; + + // Derived constants + const int n_receivers_; + const int n_baselines_; + + // Timings + common::NSTimer computeStatisticsWatch_; + common::NSTimer findbadAntennasWatch_; + common::NSTimer findBadStationsWatch_; + + // Data +#define STATS_TYPE float + std::vector<std::complex<STATS_TYPE>> stats_std_; + std::vector<std::complex<STATS_TYPE>> stats_sump2_; +}; + +} // namespace antennaflagger +} // namespace dp3 + +#endif // DPPP_ANTENNAFLAGGER_H diff --git a/antennaflagger/CMakeLists.txt b/antennaflagger/CMakeLists.txt new file mode 100644 index 000000000..1a000f310 --- /dev/null +++ b/antennaflagger/CMakeLists.txt @@ -0,0 +1 @@ +add_library(AntennaFlagger OBJECT AaartfaacFlagger.cc) diff --git a/antennaflagger/statistics.h b/antennaflagger/statistics.h new file mode 100644 index 000000000..be7a25d55 --- /dev/null +++ b/antennaflagger/statistics.h @@ -0,0 +1,228 @@ +#include <algorithm> +#include <vector> +#if defined(DEBUG) +#include <cassert> +#include <iostream> +#include <iomanip> +#endif + +namespace { + +template <typename T> +T compute_mean(const std::vector<T>& data, const std::vector<bool>& flags) { +#if defined(DEBUG) + assert(data.size() == flags.size()); +#endif + + T sum = 0; + size_t m = 0; + + for (size_t i = 0; i < data.size(); i++) { + if (!flags[i]) { + sum += data[i]; + m++; + } + } + + return m > 0 ? sum / m : 0; +} + +/* + * Computes the median of the data with length n, ignoring all values + * where the corresponding flag is set to True. + * This code is based on: + * https://medium.com/@nxtchg/calculating-median-without-sorting-eaa639cedb9f + * with a number of notable changes: + * - The input is an array with values of type T (instead of objects with value + * and weight) + * - Skips flagged values (NaN or Inf) + * - Handles even-length input as a special case + */ +template <typename T> +T compute_median_inplace(const std::vector<T>& data, + const std::vector<bool>& flags) { + assert(data.size() == flags.size()); + + T below = 0; + T above = 0; + T equal = compute_mean(data, flags); + + while (true) { + size_t n_below = 0; + size_t n_equal = 0; + size_t n_above = 0; + + for (size_t i = 0; i < data.size(); i++) { + if (flags[i]) { + continue; + } + T value = data[i]; + if (value > equal) { + if (n_above == 0 || value < above) { + above = value; + } + n_above++; + } else if (value < equal) { + if (n_below == 0 || value > below) { + below = value; + } + n_below++; + } else { + n_equal++; + } + } + + T mid = n_below + n_equal + n_above; + + if (mid > 0) { + mid /= 2; + } else { + return 0; + } + + if (mid >= n_below) { + if (mid <= n_below + n_equal) { + if (data.size() % 2 == 1) { + return equal; + } else { + if (mid == n_below) { + for (size_t j = 0; j < data.size(); j++) { + if (flags[j]) { + continue; + } + T value = data[j]; + if ((value > below) && (value <= above)) { + above = value; + } + } + return (below + above) / 2; + } else { + return (equal + above) / 2; + } + } + } + + equal = above; + } else { + equal = below; + } + } +} + +template <typename T> +T compute_median(std::vector<T> v) { + auto target = v.begin() + v.size() / 2; + std::nth_element(v.begin(), target, v.end()); + if (v.size() % 2 != 0) { + return *target; + } else { + auto before = target - 1; + std::nth_element(v.begin(), before, v.end()); + return (*target + *before) / static_cast<T>(2); + } +} + +template <typename T> +T compute_std(const std::vector<T>& data, const std::vector<bool>& flags) { +#if defined(DEBUG) + assert(data.size() == flags.size()); +#endif + + T sum = 0; + size_t m = 0; + + T mean = compute_mean(data, flags); + + for (size_t i = 0; i < data.size(); i++) { + if (!flags[i]) { + T value = mean - data[i]; + sum += value * value; + m++; + } + } + + return m > 0 ? std::sqrt(sum / m) : 0; +} + +template <typename T> +T compute_sump2(const std::vector<T>& data, const std::vector<bool>& flags) { + assert(data.size() == flags.size()); + + T sum = 0; + + for (size_t i = 0; i < data.size(); i++) { + if (!flags[i]) { + T value = data[i]; + sum += value * value; + } + } + + return sum; +} + +template <typename T> +inline bool is_valid(T x) { + bool is_nan = x != x; + bool is_finite = x != 2 * x; + return !is_nan && is_finite; +} + +template <typename T> +inline bool is_valid(const std::complex<T>& x) { + return is_valid(x.real()) && is_valid(x.imag()); +} + +template <typename T> +void flag_invalid(const std::vector<T>& data, std::vector<bool>& flags) { +#if defined(DEBUG) + assert(data.size() == flags.size()); +#endif + + for (size_t i = 0; i < data.size(); i++) { + if (!is_valid(data[i])) { + flags[i] = true; + } + } +} + +template <typename T> +std::vector<bool> find_outliers(float sigma, int maxiters, + const std::vector<T>& data) { + std::vector<bool> flags(data.size(), {false}); + flag_invalid(data, flags); + + for (int i = 0; i < maxiters; i++) { + // Compute median + T median = compute_median(data); + + // Compute standard deviation + T std = compute_std(data, flags); + + // Find outliers + size_t m = 0; + T lower_bound = median - (sigma * std); + T upper_bound = median + (sigma * std); + + for (size_t j = 0; j < data.size(); j++) { + T value = data[j]; + if (!flags[j] && (value < lower_bound || value > upper_bound)) { + flags[j] = true; + m++; +#if defined(DEBUG) + std::cout << "flagging data[" << std::setw(2) << std::scientific << j + << "] = " << value << " (bounds: [" << lower_bound << "," + << upper_bound << "]" << std::endl; +#endif + } + } + + // Exit early when no new outliers were found + if (m == 0) { + break; + } + } + + return flags; +} + +} // anonymous namespace \ No newline at end of file diff --git a/common/baseline_utils.h b/common/baseline_utils.h index 70f320b4b..a560c4562 100644 --- a/common/baseline_utils.h +++ b/common/baseline_utils.h @@ -6,17 +6,69 @@ #ifndef LOFAR_COMMON_BASELINE_UTILS_H_ #define LOFAR_COMMON_BASELINE_UTILS_H_ +#include <algorithm> +#include <cmath> #include <cstdint> #include <utility> namespace dp3 { namespace common { +/* + * Baselines are typically ordered in either column-major, or + * in row-major order. In both cases, for any baseline + * (antenna1,antenna2): antenna1 < antenna2 (i.e. there is a + * baseline (0,1), but not its mirror (1,0)). + * + * This is an example of column-major ordering with + * 5 antennas in the format "index:(antenna1,antenna2)": + * 0:(0,0) + * 1:(0,1) 2:(1,1) + * 3:(0,2) 4:(1,2) 5:(2,2) + * 6:(0,3) 7:(1,3) 8:(2,3) 9:(3,3) + * 10:(0,4) 11:(1,4) 12:(2,4) 13:(3,4) 14:(4,4) + * The first column has all the baselines with antenna1 == 0 + * + * And similarly for the row-major ordering: + * 0:(0,0) 1:(0,1) 2:(0,2) 3:(0,3) 4:(0,4) + * 5:(1,1) 6:(1,2) 7:(1,3) 8:(1,4) + * 9:(2,2) 10:(2,3) 11:(2,4) + * 12:(3,3) 13:(3,4) + * 14:(4,4) + * The first row has all the baselines with antenna1 == 0 + * + * Next, we list all baselines that contain a given antenna. + * The format is as follows: "antenna: [indices] [baselines]". + * + * First, column-major ordering: + * 0: [ 0, 1, 3, 6, 10] [(0,0), (0,1), (0,2), (0,3), (0,4)] + * 1: [ 1, 2, 4, 7, 11] [(0,1), (1,1), (1,2), (1,3), (1,4)] + * 2: [ 3, 4, 5, 8, 12] [(0,2), (1,2), (2,2), (2,3), (2,4)] + * 3: [ 6, 7, 8, 9, 13] [(0,3), (1,3), (2,3), (3,3), (3,4)] + * 4: [10, 11, 12, 13, 14] [(0,4), (1,4), (2,4), (3,4), (4,4)] + * + * Next, row-mjaor ordering: + * 0: [0, 1, 2, 3, 4] [(0,0], (0,1), (0,2), (0,3), (0,4)] + * 1: [1, 5, 6, 7, 8] [(0,1], (1,1), (1,2), (1,3), (1,4)] + * 2: [2, 6, 9, 10, 11] [(0,2], (1,2), (2,2), (2,3), (2,4)] + * 3: [3, 7, 10, 12, 13] [(0,3], (1,3), (2,3), (3,3), (3,4)] + * 4: [4, 8, 11, 13, 14] [(0,4], (1,4), (2,4), (3,4), (4,4)] + * + * This file contains several helper functions to go from index + * to baseline, from baseline to index or to get all the baselines + * for a given antenna. + */ + +enum BaselineOrder { COLUMN_MAJOR, ROW_MAJOR }; + /* * Computes the number of baselines, including auto-correlations. */ inline int ComputeNBaselines(int nAntennas) { - return nAntennas * (nAntennas - 1) / 2 + nAntennas; + int nCorrelations = (nAntennas * (nAntennas - 1)) / 2; + int nAutocorrelations = nAntennas; + int nBaselines = nCorrelations + nAutocorrelations; + return nBaselines; } /* @@ -24,50 +76,67 @@ inline int ComputeNBaselines(int nAntennas) { * number of antennas. It assumes that auto-correlations of antennas are * included. */ -inline int ComputeBaselineIndex(int a, int b, int nAntennas) { - return (a * nAntennas) + b - a * (a + 1) / 2; +inline int ComputeBaselineIndex(int a, int b, int nAntennas, + BaselineOrder order) { + int row = std::min(a, b); + int col = std::max(a, b); + if (order == ROW_MAJOR) { + return int((row * nAntennas) + col - row * (row + 1) / 2); + } else { + return int((col * (col + 1)) / 2) + row; + } } /* * Computes the antenna pair (baseline) given the baseline index. It assumes * that auto-correlations of antennas are included. */ -inline std::pair<int, int> ComputeBaseline(int baselineIndex, int nAntennas) { - const int nBaselines = ComputeNBaselines(nAntennas); - // Compute first antenna for baselineIndex - const int antenna1 = - nAntennas - (-1 + std::sqrt(8 * (nBaselines - baselineIndex))) / 2; +inline std::pair<int, int> ComputeBaseline(int baselineIndex, int nAntennas, + BaselineOrder order) { + if (order == ROW_MAJOR) { + // Compute number of baselines + const int nBaselines = ComputeNBaselines(nAntennas); + + // Compute antenna + const int antenna = + (1 + std::sqrt(1 + 8 * (nBaselines - baselineIndex - 1))) / 2; - // Compute the number of baselines (antenna1,*) - const int nBaselinesAntenna1 = - ((nAntennas * 2 + 1) * antenna1 - (antenna1 * antenna1)) / 2; + // Compute antenna1 + const int antenna1 = nAntennas - antenna; - // Compute the second antenna for baselineIndex - const int antenna2 = baselineIndex - nBaselinesAntenna1 + antenna1; + // Compute number of baselines (a,b) with a < antenna + const int n = nBaselines - nAntennas - (antenna * (antenna - 1)) / 2; - return {antenna1, antenna2}; + // Compute antenna2 + const int antenna2 = baselineIndex - n; + + // Return baseline + return {antenna1, antenna2}; + } else { + // Compute antenna2 + const int antenna2 = (1 + std::sqrt(1 + 8 * baselineIndex)) / 2; + + // Compute the number of baselines (a,b) with b < antenna2 + const int n = int((antenna2 * (antenna2 - 1)) / 2); + + // Compute antenna1 + const int antenna1 = baselineIndex - n; + + // Return baseline + return {antenna1, antenna2 - 1}; + } } /* * Computes the list of baselines for a given antenna. It assumes that * auto-correlations are included. */ -inline std::vector<int> ComputeBaselineList(int antenna, int nAntennas) { +inline std::vector<int> ComputeBaselineList(int antenna, int nAntennas, + BaselineOrder order) { std::vector<int> indices(nAntennas); - for (int i = antenna; i < nAntennas; i++) { - // Compute the number of baselines up to antenna i - int m = int((i * (i + 1)) / 2); - - if (i == antenna) { - // Add all baselines (j, antenna) - for (int j = 0; j < antenna + 1; j++) { - indices[j] = m + j; - } - } else { - // Add baseline (antenna, i) - indices[i] = m + antenna; - } + for (int i = 0; i < nAntennas; i++) { + indices[i] = ComputeBaselineIndex(antenna, i, nAntennas, order); } return indices; diff --git a/common/test/unit/tBaselineUtils.cc b/common/test/unit/tBaselineUtils.cc new file mode 100644 index 000000000..a9edc515d --- /dev/null +++ b/common/test/unit/tBaselineUtils.cc @@ -0,0 +1,118 @@ +#include "../../baseline_utils.h" + +#include <boost/test/unit_test.hpp> + +BOOST_AUTO_TEST_SUITE(baselineutils) + +BOOST_AUTO_TEST_CASE(nbaselines) { + for (int nAntennas = 2; nAntennas < 10; nAntennas++) { + int nBaselines = nAntennas * (nAntennas - 1) / 2 + nAntennas; + BOOST_TEST(nBaselines == dp3::common::ComputeNBaselines(nAntennas)); + } +} + +/* + * Helper function to test the row-major baseline utilities + */ +void testComputeBaselineRowMajor(int nAntennas) { + const dp3::common::BaselineOrder order = + dp3::common::BaselineOrder::ROW_MAJOR; + + // Iterate all baselines + int i = 0; + for (int antenna1 = 0; antenna1 < nAntennas; antenna1++) { + for (int antenna2 = antenna1; antenna2 < nAntennas; antenna2++) { + // Check reference and computed baseline + std::pair<int, int> baseline_reference = {antenna1, antenna2}; + std::pair<int, int> baseline = + dp3::common::ComputeBaseline(i, nAntennas, order); + BOOST_TEST(baseline_reference.first == baseline.first); + BOOST_TEST(baseline_reference.second == baseline.second); + + // Check baseline indices + int baselineIndex = dp3::common::ComputeBaselineIndex(antenna1, antenna2, + nAntennas, order); + BOOST_TEST(i == baselineIndex); + i++; + } + + // Check baseline indices for antenna1, it should be nAntennas + std::vector<int> baselineIndices = + dp3::common::ComputeBaselineList(antenna1, nAntennas, order); + BOOST_TEST(int(baselineIndices.size()) == nAntennas); + + // Use the baseline indices to lookup the baseline and check whether + // antenna1 occurs exactly nAntennas + 1 times (correlations + 1 // + // autocorrelation). + int count = 0; + for (int antenna2 = 0; antenna2 < nAntennas; antenna2++) { + const int baselineIndex = dp3::common::ComputeBaselineIndex( + antenna1, antenna2, nAntennas, order); + BOOST_TEST(baselineIndices[antenna2] == baselineIndex); + std::pair<int, int> baseline = + dp3::common::ComputeBaseline(baselineIndex, nAntennas, order); + count += baseline.first == antenna1; + count += baseline.second == antenna1; + } + BOOST_TEST(count == nAntennas + 1); + } +} + +BOOST_AUTO_TEST_CASE(rowmajor) { + testComputeBaselineRowMajor(5); + testComputeBaselineRowMajor(10); +} + +/* + * Helper function to test the column-major baseline utilities + */ +void testComputeBaselineColumnMajor(int nAntennas) { + const dp3::common::BaselineOrder order = + dp3::common::BaselineOrder::COLUMN_MAJOR; + + // Iterate all baselines + int i = 0; + for (int antenna2 = 0; antenna2 < nAntennas; antenna2++) { + for (int antenna1 = 0; antenna1 <= antenna2; antenna1++) { + // Check reference and computed baseline + std::pair<int, int> baseline_reference = {antenna1, antenna2}; + std::pair<int, int> baseline = + dp3::common::ComputeBaseline(i, nAntennas, order); + BOOST_TEST(baseline_reference.first == baseline.first); + BOOST_TEST(baseline_reference.second == baseline.second); + + // Check baseline indices + int baselineIndex = dp3::common::ComputeBaselineIndex(antenna1, antenna2, + nAntennas, order); + BOOST_TEST(i == baselineIndex); + i++; + } + + // Check baseline indices for antenna2, it should be nAntennas + std::vector<int> baselineIndices = + dp3::common::ComputeBaselineList(antenna2, nAntennas, order); + BOOST_TEST(int(baselineIndices.size()) == nAntennas); + + // Use the baseline indices to lookup the baseline and check whether + // antenna1 occurs exactly nAntennas + 1 times (correlations + 1 // + // autocorrelation). + int count = 0; + for (int antenna1 = 0; antenna1 < nAntennas; antenna1++) { + const int baselineIndex = dp3::common::ComputeBaselineIndex( + antenna1, antenna2, nAntennas, order); + BOOST_TEST(baselineIndices[antenna1] == baselineIndex); + std::pair<int, int> baseline = + dp3::common::ComputeBaseline(baselineIndex, nAntennas, order); + count += baseline.first == antenna2; + count += baseline.second == antenna2; + } + BOOST_TEST(count == nAntennas + 1); + } +} + +BOOST_AUTO_TEST_CASE(colmajor) { + testComputeBaselineColumnMajor(5); + testComputeBaselineColumnMajor(10); +} + +BOOST_AUTO_TEST_SUITE_END() \ No newline at end of file diff --git a/docs/schemas/AntennaFlagger.yml b/docs/schemas/AntennaFlagger.yml new file mode 100644 index 000000000..130c95e87 --- /dev/null +++ b/docs/schemas/AntennaFlagger.yml @@ -0,0 +1,24 @@ +description: >- + Flag all baselines containing bad antennas. The antennas can be specified by their number (see `selection`), and/or can be determined automatically based on (per time slot) statistics when `flagging=True`. +inputs: + step_name: + type: string + doc: unique name for the step `.` + default: antennaflagger + type: + type: string + doc: Case-insensitive step type; must be 'antennaflagger' `.` + default: antennaflagger + selection: + default: "[]" + type: string? + doc: >- + The antenna numbers (not names) to flag, currently the following selection are implemented: + Single antennas (e.g. 0,3,15,4), antenna ranges (e.g. 1~15), + autocorrelations for a selection (e.g. 0~3&&& flags all autocorrelations of antennas 0, 1, 2, 3). + Multiple selections can be separated with a semicolon (;) `.` + detect_outliers: + type: boolean + doc: >- + If true, outlier detection on antenna and station level is enabled and affected baselines are flagged. + Flagging is applied to all antennas, regardless of whether they are in `selection` or not `.` diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index 0011914be..a9c8a1287 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -1,10 +1,20 @@ +// AntennaFlagger.cc: DPPP step class to flag data with a baseline selection +// and/or using outlier detection on station and antenna level. +// +// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy) +// SPDX-License-Identifier: GPL-3.0-or-later #include "AntennaFlagger.h" -#include "../common/ParameterSet.h" -#include "../common/baseline_utils.h" -#include <boost/algorithm/string.hpp> #include <regex> +#include <set> +#include <string> +#include <utility> +#include <vector> + +#include <boost/algorithm/string.hpp> + +#include "../common/ParameterSet.h" namespace { /** @@ -56,6 +66,11 @@ std::vector<std::pair<int, int>> convertSelection( // Convert all antenna numbers from input string to a vector of ints. std::vector<int> antennas; for (const std::string& antenna : antenna_strings) { + // Skip empty inputs + if (antenna == "") { + continue; + } + // This regex searches for a range, low~high (e.g. 0~20). std::regex range("(\\d+)~(\\d+)"); std::smatch matches; @@ -99,20 +114,21 @@ std::vector<std::pair<int, int>> convertSelection( } casacore::Matrix<bool> convert(const std::string& selectionString, - uint nAntennas) { + size_t nAntennas) { assert(nAntennas != 0); casacore::Matrix<bool> bl(nAntennas, nAntennas, false); std::vector<std::pair<int, int>> selection = - convertSelection(selectionString, int(nAntennas)); + convertSelection(selectionString, static_cast<int>(nAntennas)); for (std::pair<int, int> p : selection) { bl(p.first, p.second) = true; bl(p.second, p.first) = true; } + return bl; } -} // namespace +} // anonymous namespace namespace dp3 { namespace steps { @@ -126,22 +142,47 @@ void AntennaFlagger::show(std::ostream& ostream) const { AntennaFlagger::AntennaFlagger(InputStep* input, const common::ParameterSet& parset, - const std::string& prefix) - : name_(prefix) { + const string& prefix) { initializationTimer_.start(); - // Parse arguments std::string selectionString = parset.getString(prefix + "selection", string()); - - if (selectionString.empty()) { - throw std::runtime_error("Selection needs to be set for AntennaFlagger."); - } + this->doDetectOutliers_ = parset.getBool(prefix + "detect_outliers", false); // Parse our selection string into a matrix of booleans. - n_antennas_ = input->getInfo().nantenna(); - selectionString_ = selectionString; - selection_ = convert(selectionString_, n_antennas_); + this->n_antennas_ = input->getInfo().nantenna(); + this->n_correlations_ = input->getInfo().ncorr(); + this->n_channels_ = input->getInfo().nchan(); + this->selectionString_ = selectionString; + this->selection_ = convert(selectionString, n_antennas_); + this->name_ = prefix; + + // Initialize flagger class + if (doDetectOutliers_) { + // Find out if the observation has LBA, HBA or AARTFAAC-12 data. + std::string antennaSet(input->getInfo().antennaSet()); + std::vector<std::string> antennaNames(input->getInfo().antennaNames()); + size_t n_receivers_per_station = 0; + size_t n_stations = 0; + if (antennaSet.substr(0, 3) == "LBA") { + // LOFAR LBA + n_receivers_per_station = 1; + n_stations = input->getInfo().nantenna(); + } else if (antennaSet.substr(0, 3) == "HBA") { + // LOFAR HBA + n_receivers_per_station = 1; + n_stations = input->getInfo().nantenna(); + } else if (antennaNames[0].substr(0, 3) == "A12") { + // AARTFAAC-12 + n_receivers_per_station = 48; + n_stations = 12; + } else { + throw std::runtime_error("Could not determine observation type."); + } + + this->flagger_.reset(new dp3::antennaflagger::AntennaFlagger( + n_stations, n_receivers_per_station, n_channels_, n_correlations_)); + } initializationTimer_.stop(); } @@ -151,22 +192,62 @@ bool AntennaFlagger::process(const base::DPBuffer& inBuffer) { buffer_.copy(inBuffer); - // Flags are in format: corr,chan,baseline. - casacore::Cube<bool>& flags = buffer_.getFlags(); + // Flags are in format: corr,chan,baseline. + auto flags = buffer_.getFlags(); + + /* + * Do additional flagging on the input data based on outlier statistics. + */ + if (doDetectOutliers_) { + computationTimer_.stop(); + flaggingTimer_.start(); + std::vector<std::complex<float>> data; + data.assign(buffer_.getData().begin(), buffer_.getData().end()); + flagger_->computeStats(data); + + // Find outliers both at antenna and station level + std::vector<int> flaggedAntennas1 = flagger_->findBadAntennas(); + std::vector<int> flaggedAntennas2 = flagger_->findBadStations(); + + // Combine vectors into a single set. + std::set<int> flaggedAntennas(flaggedAntennas2.begin(), + flaggedAntennas2.end()); + flaggedAntennas.insert(flaggedAntennas1.begin(), flaggedAntennas1.end()); + +#if defined(DEBUG) + // Output a list of the flagged antennas + for (int ant : collection) { + std::cout << ant << ","; + } + std::cout << std::endl; +#endif + + // Flag all antennas + for (int antenna : flaggedAntennas) { + for (size_t i = 0; i < n_antennas_; i++) { + this->selection_(antenna, i) = true; + this->selection_(i, antenna) = true; + } + } + + flaggingTimer_.stop(); + computationTimer_.start(); + } for (size_t correlation = 0; correlation < flags.nrow(); correlation++) { for (size_t channel = 0; channel < flags.ncolumn(); channel++) { for (size_t baseline = 0; baseline < flags.nplane(); baseline++) { - const std::pair<size_t, size_t> antennas = - common::ComputeBaseline(baseline, n_antennas_); + const std::pair<size_t, size_t> antennas = dp3::common::ComputeBaseline( + baseline, n_antennas_, dp3::common::BaselineOrder::ROW_MAJOR); flags(correlation, channel, baseline) = selection_(antennas.first, antennas.second); } } } - computationTimer_.stop(); + buffer_.setFlags(flags); + computationTimer_.stop(); getNextStep()->process(buffer_); return true; @@ -174,16 +255,20 @@ bool AntennaFlagger::process(const base::DPBuffer& inBuffer) { void AntennaFlagger::showTimings(std::ostream& ostream, double duration) const { double initTime = initializationTimer_.getElapsed(); + double flagTime = flaggingTimer_.getElapsed(); double computeTime = computationTimer_.getElapsed(); - double totTime = initTime + computeTime; + double totTime = initTime + flagTime + computeTime; ostream << " "; base::FlagCounter::showPerc1(ostream, totTime, duration); ostream << " AntennaFlagger " << name_ << "\n "; base::FlagCounter::showPerc1(ostream, initTime, totTime); ostream << " of it spent in initialization.\n "; + base::FlagCounter::showPerc1(ostream, flagTime, totTime); + ostream << " of it spent in flagging.\n "; base::FlagCounter::showPerc1(ostream, computeTime, totTime); - ostream << " of it spent in computation.\n"; + ostream << " of it spent in setting flags.\n"; } + } // namespace steps -} // namespace dp3 \ No newline at end of file +} // namespace dp3 diff --git a/steps/AntennaFlagger.h b/steps/AntennaFlagger.h index b5e85ce4a..6830a1e1e 100644 --- a/steps/AntennaFlagger.h +++ b/steps/AntennaFlagger.h @@ -1,19 +1,21 @@ -#ifndef DPPP_ANTENNAFLAGGER_H -#define DPPP_ANTENNAFLAGGER_H +// AntennaFlagger.h: DPPP step class to flag data with a baseline selection +// and/or using outlier detection on station and antenna level. +// +// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy) +// SPDX-License-Identifier: GPL-3.0-or-later -#include "InputStep.h" +#ifndef DPPP_ANTENNAFLAGGERSTEP_H +#define DPPP_ANTENNAFLAGGERSTEP_H -#include "../base/DPBuffer.h" -#include "../base/BaselineSelection.h" +#include <memory> +#include <string> +#include <vector> -#include <casacore/measures/Measures/MDirection.h> +#include "../antennaflagger/AntennaFlagger.h" -namespace dp3 { -namespace common { -class ParameterSet; -class ParameterValue; -} // namespace common +#include "InputStep.h" +namespace dp3 { namespace steps { class AntennaFlagger final : public Step { public: @@ -25,16 +27,26 @@ class AntennaFlagger final : public Step { void showTimings(std::ostream& ostream, double duration) const override; private: - string name_; - base::DPBuffer buffer_; + // Data statistics size_t n_antennas_; - std::string selectionString_; + size_t n_channels_; + size_t n_correlations_; + + // Data fields + base::DPBuffer buffer_; casacore::Matrix<bool> selection_; + std::string name_; + std::string selectionString_; + bool doDetectOutliers_; + std::unique_ptr<dp3::antennaflagger::AntennaFlagger> flagger_; + + // TImers common::NSTimer initializationTimer_; common::NSTimer computationTimer_; + common::NSTimer flaggingTimer_; }; } // namespace steps } // namespace dp3 -#endif +#endif // DPPP_ANTENNAFLAGGERSTEP_H diff --git a/steps/test/unit/tAntennaFlagger.cc b/steps/test/unit/tAntennaFlagger.cc index 76caa0d6d..15b73c492 100644 --- a/steps/test/unit/tAntennaFlagger.cc +++ b/steps/test/unit/tAntennaFlagger.cc @@ -74,7 +74,7 @@ class TestInput final : public InputStep { return false; } casacore::Cube<casacore::Complex> data(n_corr__, n_chan_, n_bl_); - for (int i = 0; i < int(data.size()); ++i) { + for (size_t i = 0; i < data.size(); ++i) { data.data()[i] = casacore::Complex(i + count_ * 10, i - 10 + count_ * 6); } casacore::Matrix<double> uvw(3, n_bl_); @@ -148,12 +148,12 @@ class TestOutput : public Step { virtual void show(std::ostream&) const {} virtual void updateInfo(const DPInfo& infoIn) { info() = infoIn; - BOOST_CHECK_EQUAL(int(infoIn.origNChan()), n_chan_); - BOOST_CHECK_EQUAL(int(infoIn.nchan()), n_chan_); - BOOST_CHECK_EQUAL(int(infoIn.ntime()), n_times_); + BOOST_CHECK_EQUAL(static_cast<int>(infoIn.origNChan()), n_chan_); + BOOST_CHECK_EQUAL(static_cast<int>(infoIn.nchan()), n_chan_); + BOOST_CHECK_EQUAL(static_cast<int>(infoIn.ntime()), n_times_); BOOST_CHECK_EQUAL(infoIn.timeInterval(), 5); - BOOST_CHECK_EQUAL(int(infoIn.nchanAvg()), 1); - BOOST_CHECK_EQUAL(int(infoIn.ntimeAvg()), 1); + BOOST_CHECK_EQUAL(static_cast<int>(infoIn.nchanAvg()), 1); + BOOST_CHECK_EQUAL(static_cast<int>(infoIn.ntimeAvg()), 1); } int count_; -- GitLab From 6cb8a19d69f63ba95e5e91e4b3057fb618f2571d Mon Sep 17 00:00:00 2001 From: Tammo Jan Dijkema <dijkema@astron.nl> Date: Wed, 10 Aug 2022 12:31:13 +0000 Subject: [PATCH 03/61] Replace 'flagging' by 'detect_outliers' --- docs/schemas/AntennaFlagger.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/schemas/AntennaFlagger.yml b/docs/schemas/AntennaFlagger.yml index 130c95e87..b132e4d12 100644 --- a/docs/schemas/AntennaFlagger.yml +++ b/docs/schemas/AntennaFlagger.yml @@ -1,5 +1,5 @@ description: >- - Flag all baselines containing bad antennas. The antennas can be specified by their number (see `selection`), and/or can be determined automatically based on (per time slot) statistics when `flagging=True`. + Flag all baselines containing bad antennas. The antennas can be specified by their number (see `selection`), and/or can be determined automatically based on (per time slot) statistics when `detect_outliers=True`. inputs: step_name: type: string -- GitLab From f44fe2600f325d46fbe8aa6ea7f1eb158f25cb93 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <veenboer@astron.nl> Date: Wed, 10 Aug 2022 12:32:27 +0000 Subject: [PATCH 04/61] Revert missing empty line at end of file --- .gitmodules | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitmodules b/.gitmodules index 373086ebe..fbd993e15 100644 --- a/.gitmodules +++ b/.gitmodules @@ -9,4 +9,5 @@ url = https://github.com/pybind/pybind11.git [submodule "external/schaap-packaging"] path = external/schaap-packaging - url = https://git.astron.nl/RD/schaap-packaging.git \ No newline at end of file + url = https://git.astron.nl/RD/schaap-packaging.git + -- GitLab From 35d013f958d64249091fc104432531466156c41d Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Mon, 15 Aug 2022 09:34:37 +0200 Subject: [PATCH 05/61] Add AartfaacFlagger parameters as parset keys for AntennaFlagger --- antennaflagger/AntennaFlagger.cc | 4 ++-- antennaflagger/AntennaFlagger.h | 6 +++--- steps/AntennaFlagger.cc | 17 +++++++++++++++-- steps/AntennaFlagger.h | 7 +++++++ 4 files changed, 27 insertions(+), 7 deletions(-) diff --git a/antennaflagger/AntennaFlagger.cc b/antennaflagger/AntennaFlagger.cc index 54fd5809d..f0b81e875 100644 --- a/antennaflagger/AntennaFlagger.cc +++ b/antennaflagger/AntennaFlagger.cc @@ -310,7 +310,7 @@ void AntennaFlagger::assertStatsComputed() const { } std::vector<int> AntennaFlagger::findBadAntennas(float sigma, int maxiters, - int outlier_flag_count) { + int outlier_threshold) { assertStatsComputed(); findbadAntennasWatch_.start(); @@ -332,7 +332,7 @@ std::vector<int> AntennaFlagger::findBadAntennas(float sigma, int maxiters, // Make a list of antennas that are flagged multiple times std::vector<int> flagged_antennas2; for (int antenna = 0; antenna < n_receivers_; antenna++) { - if (flagged_count[antenna] > outlier_flag_count) { + if (flagged_count[antenna] > outlier_threshold) { flagged_antennas2.push_back(antenna); } } diff --git a/antennaflagger/AntennaFlagger.h b/antennaflagger/AntennaFlagger.h index 6860a1f1a..fcc4a00c6 100644 --- a/antennaflagger/AntennaFlagger.h +++ b/antennaflagger/AntennaFlagger.h @@ -40,9 +40,9 @@ class AntennaFlagger { void computeStats(const std::vector<std::complex<double>>& data); void computeStats(const std::vector<std::complex<float>>& data); - std::vector<int> findBadAntennas(float sigma = 3, int maxiters = 5, - int outlier_flag_count = 2); - std::vector<int> findBadStations(float sigma = 2.5, int maxiters = 5); + std::vector<int> findBadAntennas(float sigma, int maxiters, + int outlier_threshold); + std::vector<int> findBadStations(float sigma, int maxiters); void reportRuntime() const; diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index a9c8a1287..fa2b6711a 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -148,6 +148,16 @@ AntennaFlagger::AntennaFlagger(InputStep* input, std::string selectionString = parset.getString(prefix + "selection", string()); this->doDetectOutliers_ = parset.getBool(prefix + "detect_outliers", false); + this->antenna_flagging_sigma_ = + parset.getFloat(prefix + "antenna_flagging_sigma", 3); + this->antenna_flagging_maxiters_ = + parset.getInt(prefix + "antenna_flagging_maxiters", 5); + this->antenna_flagging_outlier_threshold_ = + parset.getInt(prefix + "antenna_flagging_outlier_threshold", 2); + this->station_flagging_sigma_ = + parset.getFloat(prefix + "station_flagging_sigma", 2.5); + this->station_flagging_maxiters_ = + parset.getInt(prefix + "station_flagging_maxiters", 5); // Parse our selection string into a matrix of booleans. this->n_antennas_ = input->getInfo().nantenna(); @@ -206,8 +216,11 @@ bool AntennaFlagger::process(const base::DPBuffer& inBuffer) { flagger_->computeStats(data); // Find outliers both at antenna and station level - std::vector<int> flaggedAntennas1 = flagger_->findBadAntennas(); - std::vector<int> flaggedAntennas2 = flagger_->findBadStations(); + std::vector<int> flaggedAntennas1 = flagger_->findBadAntennas( + antenna_flagging_sigma_, antenna_flagging_maxiters_, + antenna_flagging_outlier_threshold_); + std::vector<int> flaggedAntennas2 = flagger_->findBadStations( + station_flagging_sigma_, station_flagging_maxiters_); // Combine vectors into a single set. std::set<int> flaggedAntennas(flaggedAntennas2.begin(), diff --git a/steps/AntennaFlagger.h b/steps/AntennaFlagger.h index 6830a1e1e..6da2d9419 100644 --- a/steps/AntennaFlagger.h +++ b/steps/AntennaFlagger.h @@ -40,6 +40,13 @@ class AntennaFlagger final : public Step { bool doDetectOutliers_; std::unique_ptr<dp3::antennaflagger::AntennaFlagger> flagger_; + // AartfaacFlagger parameters + float antenna_flagging_sigma_; + size_t antenna_flagging_maxiters_; + size_t antenna_flagging_outlier_threshold_; + float station_flagging_sigma_; + size_t station_flagging_maxiters_; + // TImers common::NSTimer initializationTimer_; common::NSTimer computationTimer_; -- GitLab From 8b6518972cd77fc1eded3b292847efc9a23b9ffa Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Mon, 15 Aug 2022 09:48:26 +0200 Subject: [PATCH 06/61] Update AntennaFlagger documentation --- docs/schemas/AntennaFlagger.yml | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/docs/schemas/AntennaFlagger.yml b/docs/schemas/AntennaFlagger.yml index b132e4d12..5536e8451 100644 --- a/docs/schemas/AntennaFlagger.yml +++ b/docs/schemas/AntennaFlagger.yml @@ -1,5 +1,8 @@ description: >- - Flag all baselines containing bad antennas. The antennas can be specified by their number (see `selection`), and/or can be determined automatically based on (per time slot) statistics when `detect_outliers=True`. + Flag all baselines containing bad antennas. The antennas can be specified by + their number (see `selection`), and/or can be determined automatically based + on (per time slot) statistics when `detect_outliers=True`. This step is mainly + intended for processing of AARTFAAC data. inputs: step_name: type: string @@ -22,3 +25,26 @@ inputs: doc: >- If true, outlier detection on antenna and station level is enabled and affected baselines are flagged. Flagging is applied to all antennas, regardless of whether they are in `selection` or not `.` + antenna_flagging_sigma: + type: float + doc: >- + Flag an atenna when its data exceeds this threshold `.` + antenna_flagging_maxiters: + type: int + doc: >- + Maximum number of iterations to find outliers in the data `.` + antenna_flagging_outlier_threshold: + type: int + doc: >- + An antenna is flagged when the sum of squares and/or standard deviation of + its data exceeds a threshold. This parameter controls whether either one of + these thresholds or both must exceeded before the antenna is considered to + be an outlier `.` + station_flagging_sigma: + type: float + doc: >- + Flag all antennas for this station when its data exceeds this threshold `.` + station_flagging_maxiters: + type: int + doc: >- + Maximum number of iterations to find outliers in the data `.` -- GitLab From a514cafa43043a1c96df49e7bd6ba8460638aaaa Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Mon, 15 Aug 2022 10:58:34 +0200 Subject: [PATCH 07/61] Fix syntax errors in documentation --- docs/schemas/AntennaFlagger.yml | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/schemas/AntennaFlagger.yml b/docs/schemas/AntennaFlagger.yml index 5536e8451..27f82736d 100644 --- a/docs/schemas/AntennaFlagger.yml +++ b/docs/schemas/AntennaFlagger.yml @@ -28,23 +28,23 @@ inputs: antenna_flagging_sigma: type: float doc: >- - Flag an atenna when its data exceeds this threshold `.` + Flag an atenna when its data exceeds this threshold `.` antenna_flagging_maxiters: type: int doc: >- - Maximum number of iterations to find outliers in the data `.` + Maximum number of iterations to find outliers in the data `.` antenna_flagging_outlier_threshold: type: int doc: >- - An antenna is flagged when the sum of squares and/or standard deviation of - its data exceeds a threshold. This parameter controls whether either one of - these thresholds or both must exceeded before the antenna is considered to - be an outlier `.` + An antenna is flagged when the sum of squares and/or standard deviation of + its data exceeds a threshold. This parameter controls whether either one of + these thresholds or both must exceeded before the antenna is considered to + be an outlier `.` station_flagging_sigma: type: float doc: >- - Flag all antennas for this station when its data exceeds this threshold `.` + Flag all antennas for this station when its data exceeds this threshold `.` station_flagging_maxiters: type: int doc: >- - Maximum number of iterations to find outliers in the data `.` + Maximum number of iterations to find outliers in the data `.` -- GitLab From 7b6aa541632e7bc27dddfd8372c6e5d4fa9bb323 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 09:26:44 +0200 Subject: [PATCH 08/61] Rename baseline_utils.h to baseline_indices.h --- antennaflagger/AntennaFlagger.h | 2 +- common/{baseline_utils.h => baseline_indices.h} | 2 +- common/test/unit/tBaselineUtils.cc | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) rename common/{baseline_utils.h => baseline_indices.h} (99%) diff --git a/antennaflagger/AntennaFlagger.h b/antennaflagger/AntennaFlagger.h index fcc4a00c6..26dec964d 100644 --- a/antennaflagger/AntennaFlagger.h +++ b/antennaflagger/AntennaFlagger.h @@ -7,7 +7,7 @@ #include <vector> #include <complex> -#include "../common/baseline_utils.h" +#include "../common/baseline_indices.h" #include "../common/Timer.h" namespace dp3 { diff --git a/common/baseline_utils.h b/common/baseline_indices.h similarity index 99% rename from common/baseline_utils.h rename to common/baseline_indices.h index a560c4562..c959bd1c1 100644 --- a/common/baseline_utils.h +++ b/common/baseline_indices.h @@ -1,6 +1,6 @@ /* * File taken from Aartfaac2ms repository: - * https://git.astron.nl/RD/aartfaac-tools/-/blob/master/common/baseline_utils.h + * https://git.astron.nl/RD/aartfaac-tools/-/blob/master/common/baseline_indices.h */ #ifndef LOFAR_COMMON_BASELINE_UTILS_H_ diff --git a/common/test/unit/tBaselineUtils.cc b/common/test/unit/tBaselineUtils.cc index a9edc515d..c08b5e0ff 100644 --- a/common/test/unit/tBaselineUtils.cc +++ b/common/test/unit/tBaselineUtils.cc @@ -1,4 +1,4 @@ -#include "../../baseline_utils.h" +#include "../../baseline_indices.h" #include <boost/test/unit_test.hpp> -- GitLab From 9c6581c20d9a2054efeac22c3b7a31b02c48e643 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 09:27:30 +0200 Subject: [PATCH 09/61] Update include guard --- common/baseline_indices.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/common/baseline_indices.h b/common/baseline_indices.h index c959bd1c1..19392a531 100644 --- a/common/baseline_indices.h +++ b/common/baseline_indices.h @@ -3,8 +3,8 @@ * https://git.astron.nl/RD/aartfaac-tools/-/blob/master/common/baseline_indices.h */ -#ifndef LOFAR_COMMON_BASELINE_UTILS_H_ -#define LOFAR_COMMON_BASELINE_UTILS_H_ +#ifndef DP3_COMMON_BASELINE_INDICES_H_ +#define DP3_COMMON_BASELINE_INDICES_H_ #include <algorithm> #include <cmath> @@ -145,4 +145,4 @@ inline std::vector<int> ComputeBaselineList(int antenna, int nAntennas, } // namespace common } // namespace dp3 -#endif // LOFAR_COMMON_BASELINE_UTILS_H_ +#endif // DP3_COMMON_BASELINE_INDICES_H_ -- GitLab From 9a7f1fce510dde3f5302fa69e4673942476374ee Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 09:44:55 +0200 Subject: [PATCH 10/61] Update BaseLineOrder to enum class, rename members --- antennaflagger/AntennaFlagger.cc | 2 +- common/baseline_indices.h | 6 +++--- common/test/unit/tBaselineUtils.cc | 2 +- steps/AntennaFlagger.cc | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/antennaflagger/AntennaFlagger.cc b/antennaflagger/AntennaFlagger.cc index f0b81e875..9a226cacb 100644 --- a/antennaflagger/AntennaFlagger.cc +++ b/antennaflagger/AntennaFlagger.cc @@ -119,7 +119,7 @@ void compute_stats_per_antenna( for (int antenna = 0; antenna < n_receivers; antenna++) { // Find all baselines with this antenna std::vector<int> indices = dp3::common::ComputeBaselineList( - antenna, n_receivers, dp3::common::BaselineOrder::ROW_MAJOR); + antenna, n_receivers, dp3::common::BaselineOrder::kRowMajor); // Sum statistics over these baselines for (int cor = 0; cor < n_correlations; cor++) { diff --git a/common/baseline_indices.h b/common/baseline_indices.h index 19392a531..d0362fcb9 100644 --- a/common/baseline_indices.h +++ b/common/baseline_indices.h @@ -59,7 +59,7 @@ namespace common { * for a given antenna. */ -enum BaselineOrder { COLUMN_MAJOR, ROW_MAJOR }; +enum class BaselineOrder { kColumnMajor, kRowMajor }; /* * Computes the number of baselines, including auto-correlations. @@ -80,7 +80,7 @@ inline int ComputeBaselineIndex(int a, int b, int nAntennas, BaselineOrder order) { int row = std::min(a, b); int col = std::max(a, b); - if (order == ROW_MAJOR) { + if (order == BaselineOrder::kRowMajor) { return int((row * nAntennas) + col - row * (row + 1) / 2); } else { return int((col * (col + 1)) / 2) + row; @@ -93,7 +93,7 @@ inline int ComputeBaselineIndex(int a, int b, int nAntennas, */ inline std::pair<int, int> ComputeBaseline(int baselineIndex, int nAntennas, BaselineOrder order) { - if (order == ROW_MAJOR) { + if (order == BaselineOrder::kRowMajor) { // Compute number of baselines const int nBaselines = ComputeNBaselines(nAntennas); diff --git a/common/test/unit/tBaselineUtils.cc b/common/test/unit/tBaselineUtils.cc index c08b5e0ff..159f183ae 100644 --- a/common/test/unit/tBaselineUtils.cc +++ b/common/test/unit/tBaselineUtils.cc @@ -16,7 +16,7 @@ BOOST_AUTO_TEST_CASE(nbaselines) { */ void testComputeBaselineRowMajor(int nAntennas) { const dp3::common::BaselineOrder order = - dp3::common::BaselineOrder::ROW_MAJOR; + dp3::common::BaselineOrder::kRowMajorJOR; // Iterate all baselines int i = 0; diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index fa2b6711a..ea1ea2a9c 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -251,7 +251,7 @@ bool AntennaFlagger::process(const base::DPBuffer& inBuffer) { for (size_t channel = 0; channel < flags.ncolumn(); channel++) { for (size_t baseline = 0; baseline < flags.nplane(); baseline++) { const std::pair<size_t, size_t> antennas = dp3::common::ComputeBaseline( - baseline, n_antennas_, dp3::common::BaselineOrder::ROW_MAJOR); + baseline, n_antennas_, dp3::common::BaselineOrder::kRowMajor); flags(correlation, channel, baseline) = selection_(antennas.first, antennas.second); } -- GitLab From f7994d50babadf6a2aabd2340fef87d77067f86f Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 09:50:17 +0200 Subject: [PATCH 11/61] Rename variables throughout baseline_indices --- common/baseline_indices.h | 40 +++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/common/baseline_indices.h b/common/baseline_indices.h index d0362fcb9..9d884d38f 100644 --- a/common/baseline_indices.h +++ b/common/baseline_indices.h @@ -64,11 +64,11 @@ enum class BaselineOrder { kColumnMajor, kRowMajor }; /* * Computes the number of baselines, including auto-correlations. */ -inline int ComputeNBaselines(int nAntennas) { - int nCorrelations = (nAntennas * (nAntennas - 1)) / 2; - int nAutocorrelations = nAntennas; - int nBaselines = nCorrelations + nAutocorrelations; - return nBaselines; +inline int ComputeNBaselines(int n_antennas) { + int n_correlations = (n_antennas * (n_antennas - 1)) / 2; + int n_autocorrelations = n_antennas; + int n_baselines = n_correlations + n_autocorrelations; + return n_baselines; } /* @@ -76,10 +76,10 @@ inline int ComputeNBaselines(int nAntennas) { * number of antennas. It assumes that auto-correlations of antennas are * included. */ -inline int ComputeBaselineIndex(int a, int b, int nAntennas, +inline int ComputeBaselineIndex(int antenna_a, int antenna_b, int nAntennas, BaselineOrder order) { - int row = std::min(a, b); - int col = std::max(a, b); + int row = std::min(antenna_a, antenna_b); + int col = std::max(antenna_a, antenna_b); if (order == BaselineOrder::kRowMajor) { return int((row * nAntennas) + col - row * (row + 1) / 2); } else { @@ -91,36 +91,36 @@ inline int ComputeBaselineIndex(int a, int b, int nAntennas, * Computes the antenna pair (baseline) given the baseline index. It assumes * that auto-correlations of antennas are included. */ -inline std::pair<int, int> ComputeBaseline(int baselineIndex, int nAntennas, +inline std::pair<int, int> ComputeBaseline(int baseline_index, int n_antennas, BaselineOrder order) { if (order == BaselineOrder::kRowMajor) { // Compute number of baselines - const int nBaselines = ComputeNBaselines(nAntennas); + const int n_baselines = ComputeNBaselines(n_antennas); // Compute antenna const int antenna = - (1 + std::sqrt(1 + 8 * (nBaselines - baselineIndex - 1))) / 2; + (1 + std::sqrt(1 + 8 * (n_baselines - baseline_index - 1))) / 2; // Compute antenna1 - const int antenna1 = nAntennas - antenna; + const int antenna1 = n_antennas - antenna; // Compute number of baselines (a,b) with a < antenna - const int n = nBaselines - nAntennas - (antenna * (antenna - 1)) / 2; + const int n = n_baselines - n_antennas - (antenna * (antenna - 1)) / 2; // Compute antenna2 - const int antenna2 = baselineIndex - n; + const int antenna2 = baseline_index - n; // Return baseline return {antenna1, antenna2}; } else { // Compute antenna2 - const int antenna2 = (1 + std::sqrt(1 + 8 * baselineIndex)) / 2; + const int antenna2 = (1 + std::sqrt(1 + 8 * baseline_index)) / 2; // Compute the number of baselines (a,b) with b < antenna2 const int n = int((antenna2 * (antenna2 - 1)) / 2); // Compute antenna1 - const int antenna1 = baselineIndex - n; + const int antenna1 = baseline_index - n; // Return baseline return {antenna1, antenna2 - 1}; @@ -131,12 +131,12 @@ inline std::pair<int, int> ComputeBaseline(int baselineIndex, int nAntennas, * Computes the list of baselines for a given antenna. It assumes that * auto-correlations are included. */ -inline std::vector<int> ComputeBaselineList(int antenna, int nAntennas, +inline std::vector<int> ComputeBaselineList(int antenna, int n_antennas, BaselineOrder order) { - std::vector<int> indices(nAntennas); + std::vector<int> indices(n_antennas); - for (int i = 0; i < nAntennas; i++) { - indices[i] = ComputeBaselineIndex(antenna, i, nAntennas, order); + for (int i = 0; i < n_antennas; i++) { + indices[i] = ComputeBaselineIndex(antenna, i, n_antennas, order); } return indices; -- GitLab From 9d18b6fa72b77a3cfa1d1d6284cdbedc739eb35b Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 09:50:43 +0200 Subject: [PATCH 12/61] Fix typo --- common/baseline_indices.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common/baseline_indices.h b/common/baseline_indices.h index 9d884d38f..861fdc645 100644 --- a/common/baseline_indices.h +++ b/common/baseline_indices.h @@ -47,7 +47,7 @@ namespace common { * 3: [ 6, 7, 8, 9, 13] [(0,3), (1,3), (2,3), (3,3), (3,4)] * 4: [10, 11, 12, 13, 14] [(0,4), (1,4), (2,4), (3,4), (4,4)] * - * Next, row-mjaor ordering: + * Next, row-major ordering: * 0: [0, 1, 2, 3, 4] [(0,0], (0,1), (0,2), (0,3), (0,4)] * 1: [1, 5, 6, 7, 8] [(0,1], (1,1), (1,2), (1,3), (1,4)] * 2: [2, 6, 9, 10, 11] [(0,2], (1,2), (2,2), (2,3), (2,4)] -- GitLab From ba179bf10b3d10acac86c602a4dc5aee2c5c5c28 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 10:09:32 +0200 Subject: [PATCH 13/61] Use size_t for row and col --- common/baseline_indices.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/common/baseline_indices.h b/common/baseline_indices.h index 861fdc645..86960c537 100644 --- a/common/baseline_indices.h +++ b/common/baseline_indices.h @@ -78,8 +78,8 @@ inline int ComputeNBaselines(int n_antennas) { */ inline int ComputeBaselineIndex(int antenna_a, int antenna_b, int nAntennas, BaselineOrder order) { - int row = std::min(antenna_a, antenna_b); - int col = std::max(antenna_a, antenna_b); + size_t row = std::min(antenna_a, antenna_b); + size_t col = std::max(antenna_a, antenna_b); if (order == BaselineOrder::kRowMajor) { return int((row * nAntennas) + col - row * (row + 1) / 2); } else { -- GitLab From 8bb51d1e98eb377fc0edd462523676694f6a2e7b Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 10:21:32 +0200 Subject: [PATCH 14/61] Use vector reserve and push_back in ComputeBaselineList --- common/baseline_indices.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/common/baseline_indices.h b/common/baseline_indices.h index 86960c537..d2b62bbf4 100644 --- a/common/baseline_indices.h +++ b/common/baseline_indices.h @@ -10,6 +10,7 @@ #include <cmath> #include <cstdint> #include <utility> +#include <vector> namespace dp3 { namespace common { @@ -133,10 +134,11 @@ inline std::pair<int, int> ComputeBaseline(int baseline_index, int n_antennas, */ inline std::vector<int> ComputeBaselineList(int antenna, int n_antennas, BaselineOrder order) { - std::vector<int> indices(n_antennas); + std::vector<int> indices; + indices.reserve(n_antennas); for (int i = 0; i < n_antennas; i++) { - indices[i] = ComputeBaselineIndex(antenna, i, n_antennas, order); + indices.push_back(ComputeBaselineIndex(antenna, i, n_antennas, order)); } return indices; -- GitLab From f52b7df5a941e58a0fc759cdb2c8b4b95b3aa5f2 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 10:30:57 +0200 Subject: [PATCH 15/61] Remove file documentation --- steps/AntennaFlagger.cc | 3 --- 1 file changed, 3 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index ea1ea2a9c..a1dc1a859 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -1,6 +1,3 @@ -// AntennaFlagger.cc: DPPP step class to flag data with a baseline selection -// and/or using outlier detection on station and antenna level. -// // Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy) // SPDX-License-Identifier: GPL-3.0-or-later -- GitLab From 267cd8cd8f1de73c843caa873120c56a972b7cff Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 10:40:57 +0200 Subject: [PATCH 16/61] Apply Google style --- steps/AntennaFlagger.cc | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index a1dc1a859..dde7735e8 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -30,14 +30,14 @@ namespace { * @param n_antennas the total number of antennas * @return a vector of pairs in the format (ant1, ant2) */ -std::vector<std::pair<int, int>> convertSelection( - const std::string& selectionString, int nAntennas) { +std::vector<std::pair<int, int>> ConvertSelection( + const std::string& selection_string, int n_antennas) { // The list of all baselines to be flagged std::vector<std::pair<int, int>> selection; // Split selection string on antenna-sets std::vector<std::string> antenna_sets; - boost::split(antenna_sets, selectionString, boost::is_any_of(";")); + boost::split(antenna_sets, selection_string, boost::is_any_of(";")); for (std::string& antenna_set : antenna_sets) { // Split on antennas within the antenna-set. @@ -110,13 +110,13 @@ std::vector<std::pair<int, int>> convertSelection( return selection; } -casacore::Matrix<bool> convert(const std::string& selectionString, +casacore::Matrix<bool> Convert(const std::string& selectionString, size_t nAntennas) { assert(nAntennas != 0); casacore::Matrix<bool> bl(nAntennas, nAntennas, false); std::vector<std::pair<int, int>> selection = - convertSelection(selectionString, static_cast<int>(nAntennas)); + ConvertSelection(selectionString, static_cast<int>(nAntennas)); for (std::pair<int, int> p : selection) { bl(p.first, p.second) = true; @@ -161,7 +161,7 @@ AntennaFlagger::AntennaFlagger(InputStep* input, this->n_correlations_ = input->getInfo().ncorr(); this->n_channels_ = input->getInfo().nchan(); this->selectionString_ = selectionString; - this->selection_ = convert(selectionString, n_antennas_); + this->selection_ = Convert(selectionString, n_antennas_); this->name_ = prefix; // Initialize flagger class @@ -194,12 +194,12 @@ AntennaFlagger::AntennaFlagger(InputStep* input, initializationTimer_.stop(); } -bool AntennaFlagger::process(const base::DPBuffer& inBuffer) { +bool AntennaFlagger::process(const base::DPBuffer& buffer) { computationTimer_.start(); - buffer_.copy(inBuffer); + buffer_.copy(buffer); - // Flags are in format: corr,chan,baseline. + // Flags are in format: corr,chan,baseline. auto flags = buffer_.getFlags(); /* -- GitLab From 47f4a1a50d7f50b508d6db88e8e9c0e929cf88ba Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 10:43:06 +0200 Subject: [PATCH 17/61] Move low inside for loop --- steps/AntennaFlagger.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index dde7735e8..e43bf8dd9 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -75,9 +75,8 @@ std::vector<std::pair<int, int>> ConvertSelection( try { if (std::regex_search(antenna, matches, range)) { // A range like 0~20 is in the string. - int low = std::stoi(matches[1]); int high = std::stoi(matches[2]) + 1; // High is inclusive - for (; low < high; low++) { + for (int low = std::stoi(matches[1]); low < high; low++) { antennas.push_back(low); } } else { -- GitLab From 78a2e49325d088c15bbbc716415f84abd797f84a Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 10:44:33 +0200 Subject: [PATCH 18/61] Make high const --- steps/AntennaFlagger.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index e43bf8dd9..71d070282 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -75,7 +75,7 @@ std::vector<std::pair<int, int>> ConvertSelection( try { if (std::regex_search(antenna, matches, range)) { // A range like 0~20 is in the string. - int high = std::stoi(matches[2]) + 1; // High is inclusive + const int high = std::stoi(matches[2]) + 1; // High is inclusive for (int low = std::stoi(matches[1]); low < high; low++) { antennas.push_back(low); } -- GitLab From 52e48635729aab0357b1618626861df77c8a3607 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 10:46:13 +0200 Subject: [PATCH 19/61] Replace std::endl by '\n' --- antennaflagger/AntennaFlagger.cc | 4 ++-- steps/AntennaFlagger.cc | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/antennaflagger/AntennaFlagger.cc b/antennaflagger/AntennaFlagger.cc index 9a226cacb..41b8a2451 100644 --- a/antennaflagger/AntennaFlagger.cc +++ b/antennaflagger/AntennaFlagger.cc @@ -191,7 +191,7 @@ void find_bad_antennas_(int n_stations, int n_receivers_per_station, if (n_outliers > 0) { std::cout << "found " << n_outliers << " outlier antenna(s) for correlation = " << cor - << ", station = " << station << std::endl; + << ", station = " << station << '\n'; } #endif @@ -275,7 +275,7 @@ void AntennaFlagger::reportRuntime() const { std::cout << "statistics: " << computeStatisticsWatch_.getElapsed() << ", bad_antennas: " << findbadAntennasWatch_.getElapsed() << ", bad_stations: " << findBadStationsWatch_.getElapsed() - << std::endl; + << '\n'; } void AntennaFlagger::computeStats( diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index 71d070282..db1b4e67b 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -83,7 +83,7 @@ std::vector<std::pair<int, int>> ConvertSelection( antennas.push_back(std::stoi(antenna)); } } catch (std::invalid_argument& e) { - std::cerr << "Could not parse: " << antenna << std::endl; + std::cerr << "Could not parse: " << antenna << '\n'; continue; } } @@ -228,7 +228,7 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { for (int ant : collection) { std::cout << ant << ","; } - std::cout << std::endl; + std::cout << '\n'; #endif // Flag all antennas -- GitLab From 18414b4c24e4447b9c3e51905eacee0b0e2d6bea Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 10:47:27 +0200 Subject: [PATCH 20/61] Remove redundant comments --- steps/AntennaFlagger.cc | 3 --- 1 file changed, 3 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index db1b4e67b..304644dee 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -88,15 +88,12 @@ std::vector<std::pair<int, int>> ConvertSelection( } } - // Iterate antennas and create tuple baselines and insert into vector. if (auto_correlation) { - // Ignore autocorrelated antenna indices. for (int antenna : antennas) { selection.emplace_back(antenna, antenna); } } - // Ignore both baseline combinations for the antennas in the vector. if (cross_correlation) { for (int antenna1 : antennas) { for (int antenna2 = 0; antenna2 < antenna1; antenna2++) { -- GitLab From 76152d6126ef689a4925c30f431ac743b9ecaaa0 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 10:53:04 +0200 Subject: [PATCH 21/61] Use vector reserve for selection --- steps/AntennaFlagger.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index 304644dee..c20057729 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -89,6 +89,7 @@ std::vector<std::pair<int, int>> ConvertSelection( } if (auto_correlation) { + selection.reserve(antennas.size()); for (int antenna : antennas) { selection.emplace_back(antenna, antenna); } @@ -96,6 +97,7 @@ std::vector<std::pair<int, int>> ConvertSelection( if (cross_correlation) { for (int antenna1 : antennas) { + selection.reserve(selection.size() + antenna1); for (int antenna2 = 0; antenna2 < antenna1; antenna2++) { selection.emplace_back(antenna1, antenna2); } -- GitLab From aaa5247423b30901e5e356633c7b8000dadc49da Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 10:59:30 +0200 Subject: [PATCH 22/61] Use preicrement --- antennaflagger/AntennaFlagger.cc | 47 ++++++++++++++++---------------- 1 file changed, 23 insertions(+), 24 deletions(-) diff --git a/antennaflagger/AntennaFlagger.cc b/antennaflagger/AntennaFlagger.cc index 41b8a2451..09a9e4ab7 100644 --- a/antennaflagger/AntennaFlagger.cc +++ b/antennaflagger/AntennaFlagger.cc @@ -44,14 +44,14 @@ void compute_stats_per_baseline( stats_std->resize(n_out); stats_sump2->resize(n_out); - for (size_t bl = 0; bl < n_baselines; bl++) { - for (size_t cor = 0; cor < n_correlations; cor++) { + for (size_t bl = 0; bl < n_baselines; ++bl) { + for (size_t cor = 0; cor < n_correlations; ++cor) { // Compute sum, sump2 OutType sum_real = 0; OutType sum_imag = 0; OutType sump2_real = 0; OutType sump2_imag = 0; - for (size_t chan = 0; chan < n_channels; chan++) { + for (size_t chan = 0; chan < n_channels; ++chan) { size_t idx = data_index(n_channels, n_correlations, bl, chan, cor); std::complex<InType> value = data[idx]; sum_real += value.real(); @@ -67,7 +67,7 @@ void compute_stats_per_baseline( // Compute standard deviation OutType sum2_real = 0; OutType sum2_imag = 0; - for (size_t chan = 0; chan < n_channels; chan++) { + for (size_t chan = 0; chan < n_channels; ++chan) { size_t idx = data_index(n_channels, n_correlations, bl, chan, cor); std::complex<InType> value = data[idx]; OutType diff_real = mean_real - value.real(); @@ -116,13 +116,13 @@ void compute_stats_per_antenna( stats_std_grouped->resize(n_receivers * n_correlations); stats_sump2_grouped->resize(n_receivers * n_correlations); - for (int antenna = 0; antenna < n_receivers; antenna++) { + for (int antenna = 0; antenna < n_receivers; ++antenna) { // Find all baselines with this antenna std::vector<int> indices = dp3::common::ComputeBaselineList( antenna, n_receivers, dp3::common::BaselineOrder::kRowMajor); // Sum statistics over these baselines - for (int cor = 0; cor < n_correlations; cor++) { + for (int cor = 0; cor < n_correlations; ++cor) { for (const int& bl : indices) { size_t idx_src = cor * n_baselines + bl; size_t idx_dst = cor * n_receivers + antenna; @@ -155,17 +155,17 @@ void find_bad_antennas_(int n_stations, int n_receivers_per_station, std::vector<std::vector<int>> bad_stations(n_stations); // Check XY, YX and YY - for (int cor = 1; cor < n_correlations; cor++) { + for (int cor = 1; cor < n_correlations; ++cor) { // The check sum is the list of statistics for the current correlation // - in: T2 check_sum[n_correlations][n_receivers] const std::complex<T>* check_sum = &stats[cor * n_receivers]; - for (int station = 0; station < n_stations; station++) { + for (int station = 0; station < n_stations; ++station) { // Make a list of all antennas for the current station, // and select the corresponding values from check_sum. std::vector<int> station_indices(n_receivers_per_station); std::vector<T> check_sum_station(n_receivers_per_station * 2); - for (int i = 0; i < n_receivers_per_station; i++) { + for (int i = 0; i < n_receivers_per_station; ++i) { int station_index = station * n_receivers_per_station + i; station_indices[i] = station_index; T real = std::log(check_sum[station_index].real()); @@ -184,7 +184,7 @@ void find_bad_antennas_(int n_stations, int n_receivers_per_station, // Count number of outliers int n_outliers = 0; - for (int i = 0; i < n_flags; i++) { + for (int i = 0; i < n_flags; ++i) { n_outliers += flags[i]; } #if defined(DEBUG) @@ -199,7 +199,7 @@ void find_bad_antennas_(int n_stations, int n_receivers_per_station, const int n = bad_stations[station].size(); bad_stations[station].resize(n + n_outliers); int j = 0; - for (int i = 0; i < n_receivers_per_station * 2; i++) { + for (int i = 0; i < n_receivers_per_station * 2; ++i) { if (flags[i]) { int station_index = station * n_receivers_per_station + (i % n_receivers_per_station); @@ -210,7 +210,7 @@ void find_bad_antennas_(int n_stations, int n_receivers_per_station, } // Add the bad stations to the global list of bad stations - for (int station = 0; station < n_stations; station++) { + for (int station = 0; station < n_stations; ++station) { all_bad_stations->insert(all_bad_stations->end(), bad_stations[station].begin(), bad_stations[station].end()); @@ -230,7 +230,7 @@ std::vector<std::complex<T>> compute_means(int n_stations, // Split the real and imaginary parts of the check sum std::vector<T> check_sum_real(n_receivers); std::vector<T> check_sum_imag(n_receivers); - for (int antenna = 0; antenna < n_receivers; antenna++) { + for (int antenna = 0; antenna < n_receivers; ++antenna) { check_sum_real[antenna] = check_sum[antenna].real(); check_sum_imag[antenna] = check_sum[antenna].imag(); } @@ -239,12 +239,12 @@ std::vector<std::complex<T>> compute_means(int n_stations, T median_real = compute_median(check_sum_real); T median_imag = compute_median(check_sum_imag); - for (int station = 0; station < n_stations; station++) { + for (int station = 0; station < n_stations; ++station) { // Make a list of all antennas for the current station, // and select the corresponding values from check_sum std::vector<T> check_sum_station_real(n_receivers_per_station); std::vector<T> check_sum_station_imag(n_receivers_per_station); - for (int i = 0; i < n_receivers_per_station; i++) { + for (int i = 0; i < n_receivers_per_station; ++i) { int station_index = station * n_receivers_per_station + i; check_sum_station_real[i] = check_sum_real[station_index]; check_sum_station_imag[i] = check_sum_imag[station_index]; @@ -274,8 +274,7 @@ namespace antennaflagger { void AntennaFlagger::reportRuntime() const { std::cout << "statistics: " << computeStatisticsWatch_.getElapsed() << ", bad_antennas: " << findbadAntennasWatch_.getElapsed() - << ", bad_stations: " << findBadStationsWatch_.getElapsed() - << '\n'; + << ", bad_stations: " << findBadStationsWatch_.getElapsed() << '\n'; } void AntennaFlagger::computeStats( @@ -331,7 +330,7 @@ std::vector<int> AntennaFlagger::findBadAntennas(float sigma, int maxiters, // Make a list of antennas that are flagged multiple times std::vector<int> flagged_antennas2; - for (int antenna = 0; antenna < n_receivers_; antenna++) { + for (int antenna = 0; antenna < n_receivers_; ++antenna) { if (flagged_count[antenna] > outlier_threshold) { flagged_antennas2.push_back(antenna); } @@ -360,10 +359,10 @@ std::vector<int> AntennaFlagger::findBadStations(float sigma, int maxiters) { means.resize(n_stations_); } - for (size_t i = 0; i < check_sums.size(); i++) { + for (size_t i = 0; i < check_sums.size(); ++i) { std::vector<std::complex<STATS_TYPE>> means = compute_means(n_stations_, n_receivers_per_station, check_sums[i]); - for (int station = 0; station < n_stations_; station++) { + for (int station = 0; station < n_stations_; ++station) { all_means[0 * check_sums.size() + i][station] = means[station].real(); all_means[1 * check_sums.size() + i][station] = means[station].imag(); } @@ -371,7 +370,7 @@ std::vector<int> AntennaFlagger::findBadStations(float sigma, int maxiters) { // Flag all stations for which all of the means are an outlier std::vector<bool> station_flags(n_stations_, {false}); - for (size_t i = 0; i < all_means.size(); i++) { + for (size_t i = 0; i < all_means.size(); ++i) { const std::vector<STATS_TYPE>& means = all_means[i]; // Find outliers @@ -380,7 +379,7 @@ std::vector<int> AntennaFlagger::findBadStations(float sigma, int maxiters) { assert(flags.size() == size_t(n_stations_)); #endif - for (int station = 0; station < n_stations_; station++) { + for (int station = 0; station < n_stations_; ++station) { station_flags[station] = (flags[station] && (i == 0 || station_flags[station])); } @@ -388,14 +387,14 @@ std::vector<int> AntennaFlagger::findBadStations(float sigma, int maxiters) { // Make a list of all the antenna indices for the flagged stations std::vector<int> flagged_antennas; - for (int station = 0; station < n_stations_; station++) { + for (int station = 0; station < n_stations_; ++station) { if (station_flags[station]) { // Make space for the new indices int n = flagged_antennas.size(); flagged_antennas.resize(n + n_receivers_per_station); // Add the new indices - for (int i = 0; i < n_receivers_per_station; i++) { + for (int i = 0; i < n_receivers_per_station; ++i) { int station_index = station * n_receivers_per_station + i; flagged_antennas[n + i] = station_index; } -- GitLab From 216e3d964e067bd1f42f5d2a689939dab8844e76 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 11:01:55 +0200 Subject: [PATCH 23/61] Add missing std namespace --- steps/AntennaFlagger.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index c20057729..64f3151d0 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -141,7 +141,7 @@ AntennaFlagger::AntennaFlagger(InputStep* input, initializationTimer_.start(); // Parse arguments std::string selectionString = - parset.getString(prefix + "selection", string()); + parset.getString(prefix + "selection", std::string()); this->doDetectOutliers_ = parset.getBool(prefix + "detect_outliers", false); this->antenna_flagging_sigma_ = parset.getFloat(prefix + "antenna_flagging_sigma", 3); -- GitLab From f534159b48eb2563fb7c4bfca9fa08e252a0699e Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 11:03:32 +0200 Subject: [PATCH 24/61] Remove this-> --- steps/AntennaFlagger.cc | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index 64f3151d0..c2584f436 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -142,25 +142,25 @@ AntennaFlagger::AntennaFlagger(InputStep* input, // Parse arguments std::string selectionString = parset.getString(prefix + "selection", std::string()); - this->doDetectOutliers_ = parset.getBool(prefix + "detect_outliers", false); - this->antenna_flagging_sigma_ = + doDetectOutliers_ = parset.getBool(prefix + "detect_outliers", false); + antenna_flagging_sigma_ = parset.getFloat(prefix + "antenna_flagging_sigma", 3); - this->antenna_flagging_maxiters_ = + antenna_flagging_maxiters_ = parset.getInt(prefix + "antenna_flagging_maxiters", 5); - this->antenna_flagging_outlier_threshold_ = + antenna_flagging_outlier_threshold_ = parset.getInt(prefix + "antenna_flagging_outlier_threshold", 2); - this->station_flagging_sigma_ = + station_flagging_sigma_ = parset.getFloat(prefix + "station_flagging_sigma", 2.5); - this->station_flagging_maxiters_ = + station_flagging_maxiters_ = parset.getInt(prefix + "station_flagging_maxiters", 5); // Parse our selection string into a matrix of booleans. - this->n_antennas_ = input->getInfo().nantenna(); - this->n_correlations_ = input->getInfo().ncorr(); - this->n_channels_ = input->getInfo().nchan(); - this->selectionString_ = selectionString; - this->selection_ = Convert(selectionString, n_antennas_); - this->name_ = prefix; + n_antennas_ = input->getInfo().nantenna(); + n_correlations_ = input->getInfo().ncorr(); + n_channels_ = input->getInfo().nchan(); + selectionString_ = selectionString; + selection_ = Convert(selectionString, n_antennas_); + name_ = prefix; // Initialize flagger class if (doDetectOutliers_) { @@ -185,7 +185,7 @@ AntennaFlagger::AntennaFlagger(InputStep* input, throw std::runtime_error("Could not determine observation type."); } - this->flagger_.reset(new dp3::antennaflagger::AntennaFlagger( + flagger_.reset(new dp3::antennaflagger::AntennaFlagger( n_stations, n_receivers_per_station, n_channels_, n_correlations_)); } @@ -233,8 +233,8 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { // Flag all antennas for (int antenna : flaggedAntennas) { for (size_t i = 0; i < n_antennas_; i++) { - this->selection_(antenna, i) = true; - this->selection_(i, antenna) = true; + selection_(antenna, i) = true; + selection_(i, antenna) = true; } } -- GitLab From ad94fa61bd9b1dc9eb3e38985e5303730c0dac54 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 11:05:07 +0200 Subject: [PATCH 25/61] Fix typo AaartfaacFlagger --- antennaflagger/AntennaFlagger.cc | 4 ++-- antennaflagger/CMakeLists.txt | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) delete mode 100644 antennaflagger/CMakeLists.txt diff --git a/antennaflagger/AntennaFlagger.cc b/antennaflagger/AntennaFlagger.cc index 09a9e4ab7..15bea4d57 100644 --- a/antennaflagger/AntennaFlagger.cc +++ b/antennaflagger/AntennaFlagger.cc @@ -20,7 +20,7 @@ inline size_t data_index(int n_channels, int n_correlations, int bl, int chan, } /* - * Helper functions for AaartfaacFlagger::find_bad_antennas + * Helper functions for AartfaacFlagger::find_bad_antennas */ template <typename InType, typename OutType> void compute_stats_per_baseline( @@ -134,7 +134,7 @@ void compute_stats_per_antenna( } // end compute_stats_per_antenna /* - * Helper functions for AaartfaacFlagger::find_bad_stations + * Helper functions for AaatfaacFlagger::find_bad_stations */ template <typename T> void find_bad_antennas_(int n_stations, int n_receivers_per_station, diff --git a/antennaflagger/CMakeLists.txt b/antennaflagger/CMakeLists.txt deleted file mode 100644 index 1a000f310..000000000 --- a/antennaflagger/CMakeLists.txt +++ /dev/null @@ -1 +0,0 @@ -add_library(AntennaFlagger OBJECT AaartfaacFlagger.cc) -- GitLab From e469c795ce18edf66b680a56e9c7559676167d81 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 11:15:39 +0200 Subject: [PATCH 26/61] Combine check for LBA and HBA --- steps/AntennaFlagger.cc | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index c2584f436..0c331ad75 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -169,12 +169,8 @@ AntennaFlagger::AntennaFlagger(InputStep* input, std::vector<std::string> antennaNames(input->getInfo().antennaNames()); size_t n_receivers_per_station = 0; size_t n_stations = 0; - if (antennaSet.substr(0, 3) == "LBA") { - // LOFAR LBA - n_receivers_per_station = 1; - n_stations = input->getInfo().nantenna(); - } else if (antennaSet.substr(0, 3) == "HBA") { - // LOFAR HBA + if (antennaSet.substr(0, 3) == "LBA" || antennaSet.substr(0, 3) == "HBA") { + // LOFAR LBA or HBA n_receivers_per_station = 1; n_stations = input->getInfo().nantenna(); } else if (antennaNames[0].substr(0, 3) == "A12") { -- GitLab From 2f69d38cdbe0462916c030d6c942d2c5573e1a64 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 11:17:36 +0200 Subject: [PATCH 27/61] Use make_unique instead of reset --- steps/AntennaFlagger.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index 0c331ad75..b3baaed0a 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -181,8 +181,8 @@ AntennaFlagger::AntennaFlagger(InputStep* input, throw std::runtime_error("Could not determine observation type."); } - flagger_.reset(new dp3::antennaflagger::AntennaFlagger( - n_stations, n_receivers_per_station, n_channels_, n_correlations_)); + flagger_ = std::make_unique<dp3::antennaflagger::AntennaFlagger>( + n_stations, n_receivers_per_station, n_channels_, n_correlations_); } initializationTimer_.stop(); -- GitLab From 381791c2beef8253ada8ee89f86723064e042727 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 11:18:56 +0200 Subject: [PATCH 28/61] Replace auto --- steps/AntennaFlagger.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index b3baaed0a..ea61e0a38 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -194,7 +194,7 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { buffer_.copy(buffer); // Flags are in format: corr,chan,baseline. - auto flags = buffer_.getFlags(); + casacore::Cube<bool> flags = buffer_.getFlags(); /* * Do additional flagging on the input data based on outlier statistics. -- GitLab From 3ac959636ad5895a936a81276f74d09b324a7e5c Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 11:19:41 +0200 Subject: [PATCH 29/61] Shorten initialization of data vector --- steps/AntennaFlagger.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index ea61e0a38..f79830cc2 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -202,8 +202,8 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { if (doDetectOutliers_) { computationTimer_.stop(); flaggingTimer_.start(); - std::vector<std::complex<float>> data; - data.assign(buffer_.getData().begin(), buffer_.getData().end()); + std::vector<std::complex<float>> data(buffer_.getData().begin(), + buffer_.getData().end()); flagger_->computeStats(data); // Find outliers both at antenna and station level -- GitLab From 2b2fe99d11c69a875ac3aaf8a466d12c77fa11d1 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 11:20:44 +0200 Subject: [PATCH 30/61] Make vectors const --- steps/AntennaFlagger.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index f79830cc2..8fc14e731 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -202,15 +202,15 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { if (doDetectOutliers_) { computationTimer_.stop(); flaggingTimer_.start(); - std::vector<std::complex<float>> data(buffer_.getData().begin(), - buffer_.getData().end()); + const std::vector<std::complex<float>> data(buffer_.getData().begin(), + buffer_.getData().end()); flagger_->computeStats(data); // Find outliers both at antenna and station level - std::vector<int> flaggedAntennas1 = flagger_->findBadAntennas( + const std::vector<int> flaggedAntennas1 = flagger_->findBadAntennas( antenna_flagging_sigma_, antenna_flagging_maxiters_, antenna_flagging_outlier_threshold_); - std::vector<int> flaggedAntennas2 = flagger_->findBadStations( + const std::vector<int> flaggedAntennas2 = flagger_->findBadStations( station_flagging_sigma_, station_flagging_maxiters_); // Combine vectors into a single set. -- GitLab From fde7936320ba64f99fd72dbec152e6f1f52e7388 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 11:21:33 +0200 Subject: [PATCH 31/61] Rename flaggedAntennas to flagged_antennas --- steps/AntennaFlagger.cc | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index 8fc14e731..d6cdb5ab0 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -207,16 +207,17 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { flagger_->computeStats(data); // Find outliers both at antenna and station level - const std::vector<int> flaggedAntennas1 = flagger_->findBadAntennas( + const std::vector<int> flagged_antennas_1 = flagger_->findBadAntennas( antenna_flagging_sigma_, antenna_flagging_maxiters_, antenna_flagging_outlier_threshold_); - const std::vector<int> flaggedAntennas2 = flagger_->findBadStations( + const std::vector<int> flagged_antennas_2 = flagger_->findBadStations( station_flagging_sigma_, station_flagging_maxiters_); // Combine vectors into a single set. - std::set<int> flaggedAntennas(flaggedAntennas2.begin(), - flaggedAntennas2.end()); - flaggedAntennas.insert(flaggedAntennas1.begin(), flaggedAntennas1.end()); + std::set<int> flagged_antennas(flagged_antennas_2.begin(), + flagged_antennas_2.end()); + flagged_antennas.insert(flagged_antennas_1.begin(), + flagged_antennas_1.end()); #if defined(DEBUG) // Output a list of the flagged antennas @@ -227,7 +228,7 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { #endif // Flag all antennas - for (int antenna : flaggedAntennas) { + for (int antenna : flagged_antennas) { for (size_t i = 0; i < n_antennas_; i++) { selection_(antenna, i) = true; selection_(i, antenna) = true; -- GitLab From fbb287ba8e9ee95a056ca3ec4b8a4b8ddfbcd75b Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 11:22:32 +0200 Subject: [PATCH 32/61] Remove file comment --- steps/AntennaFlagger.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/steps/AntennaFlagger.h b/steps/AntennaFlagger.h index 6da2d9419..0408e732d 100644 --- a/steps/AntennaFlagger.h +++ b/steps/AntennaFlagger.h @@ -1,6 +1,3 @@ -// AntennaFlagger.h: DPPP step class to flag data with a baseline selection -// and/or using outlier detection on station and antenna level. -// // Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy) // SPDX-License-Identifier: GPL-3.0-or-later -- GitLab From a6a339935e23eac351787f92f8d5729b6b6b3f55 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 11:53:51 +0200 Subject: [PATCH 33/61] Add xtl, xsimd and xtensor to cmake --- CMakeLists.txt | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index a2d539708..1e9a0ae0d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -173,6 +173,39 @@ include_directories(${Boost_INCLUDE_DIRS}) find_package(Threads REQUIRED) +include(ExternalProject) + +set(XTL_ROOT_DIR ${CMAKE_BINARY_DIR}/external/xtl) +set(XSIMD_ROOT_DIR ${CMAKE_BINARY_DIR}/external/xsimd) +set(XTENSOR_ROOT_DIR ${CMAKE_BINARY_DIR}/external/xtensor) + +ExternalProject_Add( + xtl + INSTALL_DIR ${XTL_ROOT_DIR} + GIT_REPOSITORY https://github.com/xtensor-stack/xtl.git + GIT_TAG 0.7.4 + CMAKE_ARGS -DCMAKE_INSTALL_MESSAGE=LAZY) + +ExternalProject_Add( + xsimd + DEPENDS xtl + INSTALL_DIR ${CMAKE_BINARY_DIR}/external/xsimd + GIT_REPOSITORY https://github.com/xtensor-stack/xsimd.git + GIT_TAG 8.1.0 + CMAKE_ARGS -DCMAKE_INSTALL_MESSAGE=LAZY) +ExternalProject_Add( + xtensor + DEPENDS xsimd + INSTALL_DIR ${CMAKE_BINARY_DIR}/external/xtensor + GIT_REPOSITORY https://github.com/xtensor-stack/xtensor.git + GIT_TAG 0.24.2 + CMAKE_ARGS -DXTENSOR_USE_XSIMD=ON + -DCMAKE_PREFIX_PATH=${XTL_ROOT_DIR}:${XSIMD_ROOT_DIR} + -DCMAKE_INSTALL_MESSAGE=LAZY) +set(XTL_INCLUDE_DIR ${XTL_ROOT_DIR}/include) +set(XTENSOR_INCLUDE_DIR ${XTENSOR_ROOT_DIR}/include) +include_directories(${XTL_INCLUDE_DIR} ${XTENSOR_INCLUDE_DIR}) + # === Load astron packages === find_package(AOFlagger 3.1 REQUIRED) -- GitLab From 2e9c59d24fdbab670fbebdad657eee9acc982f4c Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 15:33:15 +0200 Subject: [PATCH 34/61] Update cmake for external xtensor --- CMakeLists.txt | 89 ++++++++++++++++++++++++++++++++++---------------- 1 file changed, 60 insertions(+), 29 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1e9a0ae0d..b3f992604 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -175,35 +175,62 @@ find_package(Threads REQUIRED) include(ExternalProject) -set(XTL_ROOT_DIR ${CMAKE_BINARY_DIR}/external/xtl) -set(XSIMD_ROOT_DIR ${CMAKE_BINARY_DIR}/external/xsimd) -set(XTENSOR_ROOT_DIR ${CMAKE_BINARY_DIR}/external/xtensor) - -ExternalProject_Add( - xtl - INSTALL_DIR ${XTL_ROOT_DIR} - GIT_REPOSITORY https://github.com/xtensor-stack/xtl.git - GIT_TAG 0.7.4 - CMAKE_ARGS -DCMAKE_INSTALL_MESSAGE=LAZY) - -ExternalProject_Add( - xsimd - DEPENDS xtl - INSTALL_DIR ${CMAKE_BINARY_DIR}/external/xsimd - GIT_REPOSITORY https://github.com/xtensor-stack/xsimd.git - GIT_TAG 8.1.0 - CMAKE_ARGS -DCMAKE_INSTALL_MESSAGE=LAZY) -ExternalProject_Add( - xtensor - DEPENDS xsimd - INSTALL_DIR ${CMAKE_BINARY_DIR}/external/xtensor - GIT_REPOSITORY https://github.com/xtensor-stack/xtensor.git - GIT_TAG 0.24.2 - CMAKE_ARGS -DXTENSOR_USE_XSIMD=ON - -DCMAKE_PREFIX_PATH=${XTL_ROOT_DIR}:${XSIMD_ROOT_DIR} - -DCMAKE_INSTALL_MESSAGE=LAZY) -set(XTL_INCLUDE_DIR ${XTL_ROOT_DIR}/include) -set(XTENSOR_INCLUDE_DIR ${XTENSOR_ROOT_DIR}/include) +# Try to find xtl and xtensor from the environment first +find_path(XTL_INCLUDE_DIR xtl/xtl_config.hpp HINTS ENV ${XTL_ROOT}) +find_path(XTENSOR_INCLUDE_DIR xtensor.hpp HINTS ENV ${XTENSOR_ROOT}) + +if(${XTL_INCLUDE_DIR} STREQUAL "XTL_INCLUDE_DIR-NOTFOUND" + OR ${XTENSOR_INCLUDE_DIR} STREQUAL "XTENSOR_INCLUDE_DIR-NOTFOUND") + + # Try to find locally installed xtl and xtensor + set(XTL_ROOT_DIR ${CMAKE_BINARY_DIR}/external/xtl) + set(XSIMD_ROOT_DIR ${CMAKE_BINARY_DIR}/external/xsimd) + set(XTENSOR_ROOT_DIR ${CMAKE_BINARY_DIR}/external/xtensor) + find_path( + XTL_INCLUDE_DIR xtl/xtl_config.hpp + PATH_SUFFIXES include + HINTS ${XTL_ROOT_DIR}) + find_path( + XTENSOR_INCLUDE_DIR xtensor/xtensor.hpp + PATH_SUFFIXES include + HINTS ${XTENSOR_ROOT_DIR}) +endif() + +# If xtl and xtensor are not both found, build them locally +if(${XTL_INCLUDE_DIR} STREQUAL "XTL_INCLUDE_DIR-NOTFOUND" + OR ${XTENSOR_INCLUDE_DIR} STREQUAL "XTENSOR_INCLUDE_DIR-NOTFOUND") + set(BUILD_XTENSOR TRUE) + ExternalProject_Add( + xtl + INSTALL_DIR ${XTL_ROOT_DIR} + GIT_REPOSITORY https://github.com/xtensor-stack/xtl.git + GIT_TAG 0.7.4 + EXCLUDE_FROM_ALL TRUE + CMAKE_ARGS -DCMAKE_INSTALL_MESSAGE=LAZY) + ExternalProject_Add( + xsimd + DEPENDS xtl + INSTALL_DIR ${XSIMD_ROOT_DIR} + GIT_REPOSITORY https://github.com/xtensor-stack/xsimd.git + GIT_TAG 8.1.0 + EXCLUDE_FROM_ALL TRUE + CMAKE_ARGS -DCMAKE_INSTALL_MESSAGE=LAZY) + ExternalProject_Add( + xtensor + DEPENDS xsimd + INSTALL_DIR ${XTENSOR_ROOT_DIR} + GIT_REPOSITORY https://github.com/xtensor-stack/xtensor.git + GIT_TAG 0.24.2 + EXCLUDE_FROM_ALL TRUE + CMAKE_ARGS -DXTENSOR_USE_XSIMD=ON + -DCMAKE_PREFIX_PATH=${XTL_ROOT_DIR}:${XSIMD_ROOT_DIR} + -DCMAKE_INSTALL_MESSAGE=LAZY) + set(XTL_INCLUDE_DIR ${XTL_ROOT_DIR}/include) + set(XTENSOR_INCLUDE_DIR ${XTENSOR_ROOT_DIR}/include) +else() + # Prevent rebuild of xtensor + set(BUILD_XTENSOR FALSE) +endif() include_directories(${XTL_INCLUDE_DIR} ${XTENSOR_INCLUDE_DIR}) # === Load astron packages === @@ -460,6 +487,10 @@ add_library( steps/DemixerNew.cc steps/NullStokes.cc) +if(${BUILD_XTENSOR}) + add_dependencies(DP3_OBJ xtensor) +endif() + add_library( ParmDB OBJECT parmdb/Axis.cc -- GitLab From 0a572a84a53667733c98412c49b740a90348261f Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Thu, 18 Aug 2022 16:41:48 +0200 Subject: [PATCH 35/61] Use xtensor for selection_ member --- steps/AntennaFlagger.cc | 6 +++--- steps/AntennaFlagger.h | 4 +++- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index d6cdb5ab0..dce1892ab 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -108,10 +108,10 @@ std::vector<std::pair<int, int>> ConvertSelection( return selection; } -casacore::Matrix<bool> Convert(const std::string& selectionString, - size_t nAntennas) { +xt::xtensor<bool, 2> Convert(const std::string& selectionString, + size_t nAntennas) { assert(nAntennas != 0); - casacore::Matrix<bool> bl(nAntennas, nAntennas, false); + xt::xtensor<bool, 2> bl({nAntennas, nAntennas}, false); std::vector<std::pair<int, int>> selection = ConvertSelection(selectionString, static_cast<int>(nAntennas)); diff --git a/steps/AntennaFlagger.h b/steps/AntennaFlagger.h index 0408e732d..d309be24a 100644 --- a/steps/AntennaFlagger.h +++ b/steps/AntennaFlagger.h @@ -7,6 +7,7 @@ #include <memory> #include <string> #include <vector> +#include <xtensor.hpp> #include "../antennaflagger/AntennaFlagger.h" @@ -31,7 +32,8 @@ class AntennaFlagger final : public Step { // Data fields base::DPBuffer buffer_; - casacore::Matrix<bool> selection_; + xt::xtensor<bool, 2> selection_; + std::string name_; std::string selectionString_; bool doDetectOutliers_; -- GitLab From de1146e394524a1051f6326130aa4269423e08f9 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 08:26:39 +0200 Subject: [PATCH 36/61] Remove unneeded comments in ComputeBaseline --- common/baseline_indices.h | 20 ++------------------ 1 file changed, 2 insertions(+), 18 deletions(-) diff --git a/common/baseline_indices.h b/common/baseline_indices.h index d2b62bbf4..75f3b7664 100644 --- a/common/baseline_indices.h +++ b/common/baseline_indices.h @@ -95,35 +95,19 @@ inline int ComputeBaselineIndex(int antenna_a, int antenna_b, int nAntennas, inline std::pair<int, int> ComputeBaseline(int baseline_index, int n_antennas, BaselineOrder order) { if (order == BaselineOrder::kRowMajor) { - // Compute number of baselines const int n_baselines = ComputeNBaselines(n_antennas); - - // Compute antenna const int antenna = (1 + std::sqrt(1 + 8 * (n_baselines - baseline_index - 1))) / 2; - - // Compute antenna1 const int antenna1 = n_antennas - antenna; - - // Compute number of baselines (a,b) with a < antenna + // n is the number of baselines (a,b) with a < antenna const int n = n_baselines - n_antennas - (antenna * (antenna - 1)) / 2; - - // Compute antenna2 const int antenna2 = baseline_index - n; - - // Return baseline return {antenna1, antenna2}; } else { - // Compute antenna2 const int antenna2 = (1 + std::sqrt(1 + 8 * baseline_index)) / 2; - - // Compute the number of baselines (a,b) with b < antenna2 + // n is the number of baselines (a,b) with b < antenna2 const int n = int((antenna2 * (antenna2 - 1)) / 2); - - // Compute antenna1 const int antenna1 = baseline_index - n; - - // Return baseline return {antenna1, antenna2 - 1}; } } -- GitLab From 6b5281e203ff731172aa3d086cd998f063e564c4 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 08:40:33 +0200 Subject: [PATCH 37/61] Update comment in baseline_indices --- common/baseline_indices.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/common/baseline_indices.h b/common/baseline_indices.h index 75f3b7664..7f0a9b13f 100644 --- a/common/baseline_indices.h +++ b/common/baseline_indices.h @@ -16,10 +16,10 @@ namespace dp3 { namespace common { /* - * Baselines are typically ordered in either column-major, or - * in row-major order. In both cases, for any baseline - * (antenna1,antenna2): antenna1 < antenna2 (i.e. there is a - * baseline (0,1), but not its mirror (1,0)). + * Baselines are typically ordered in either column-major (e.g. for LOFAR), or + * in row-major order (e.g. for AARTFAAC). In both cases, for any baseline + * (antenna1,antenna2): antenna1 < antenna2 (i.e. there is a baseline (0,1), but + * not its mirror (1,0)). * * This is an example of column-major ordering with * 5 antennas in the format "index:(antenna1,antenna2)": -- GitLab From f523665f75226777239501a8468464ce3622f4e2 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 08:40:57 +0200 Subject: [PATCH 38/61] Replace int by size_t --- antennaflagger/AntennaFlagger.cc | 94 ++++++++++++++++---------------- antennaflagger/AntennaFlagger.h | 22 ++++---- common/baseline_indices.h | 47 ++++++++-------- steps/AntennaFlagger.cc | 8 +-- 4 files changed, 86 insertions(+), 85 deletions(-) diff --git a/antennaflagger/AntennaFlagger.cc b/antennaflagger/AntennaFlagger.cc index 15bea4d57..16a42c4f1 100644 --- a/antennaflagger/AntennaFlagger.cc +++ b/antennaflagger/AntennaFlagger.cc @@ -13,8 +13,8 @@ // Non-member helper functions namespace { -inline size_t data_index(int n_channels, int n_correlations, int bl, int chan, - int cor) { +inline size_t data_index(size_t n_channels, size_t n_correlations, size_t bl, + size_t chan, size_t cor) { // data[n_baselines][n_channels][n_correlations] return (bl * n_channels + chan) * n_correlations + cor; } @@ -88,8 +88,8 @@ void compute_stats_per_baseline( template <typename InType, typename OutType> void compute_stats_per_antenna( - int n_receivers, int n_baselines, int n_channels, int n_correlations, - const std::vector<std::complex<InType>>& data, + size_t n_receivers, size_t n_baselines, size_t n_channels, + size_t n_correlations, const std::vector<std::complex<InType>>& data, std::vector<std::complex<OutType>>* stats_std_grouped, std::vector<std::complex<OutType>>* stats_sump2_grouped) { // Compute stats for every baseline and correlation: @@ -116,14 +116,14 @@ void compute_stats_per_antenna( stats_std_grouped->resize(n_receivers * n_correlations); stats_sump2_grouped->resize(n_receivers * n_correlations); - for (int antenna = 0; antenna < n_receivers; ++antenna) { + for (size_t antenna = 0; antenna < n_receivers; ++antenna) { // Find all baselines with this antenna - std::vector<int> indices = dp3::common::ComputeBaselineList( + std::vector<size_t> indices = dp3::common::ComputeBaselineList( antenna, n_receivers, dp3::common::BaselineOrder::kRowMajor); // Sum statistics over these baselines - for (int cor = 0; cor < n_correlations; ++cor) { - for (const int& bl : indices) { + for (size_t cor = 0; cor < n_correlations; ++cor) { + for (const size_t& bl : indices) { size_t idx_src = cor * n_baselines + bl; size_t idx_dst = cor * n_receivers + antenna; (*stats_std_grouped)[idx_dst] += stats_std[idx_src]; @@ -137,10 +137,10 @@ void compute_stats_per_antenna( * Helper functions for AaatfaacFlagger::find_bad_stations */ template <typename T> -void find_bad_antennas_(int n_stations, int n_receivers_per_station, - int n_baselines, int n_correlations, float sigma, +void find_bad_antennas_(size_t n_stations, size_t n_receivers_per_station, + size_t n_baselines, size_t n_correlations, float sigma, int maxiters, const std::vector<std::complex<T>>& stats, - std::vector<int>* all_bad_stations) { + std::vector<size_t>* all_bad_stations) { int n_receivers = n_stations * n_receivers_per_station; // - in: stats[n_correlations][n_receivers] // - out: all_bad_stations[?] @@ -155,18 +155,18 @@ void find_bad_antennas_(int n_stations, int n_receivers_per_station, std::vector<std::vector<int>> bad_stations(n_stations); // Check XY, YX and YY - for (int cor = 1; cor < n_correlations; ++cor) { + for (size_t cor = 1; cor < n_correlations; ++cor) { // The check sum is the list of statistics for the current correlation // - in: T2 check_sum[n_correlations][n_receivers] const std::complex<T>* check_sum = &stats[cor * n_receivers]; - for (int station = 0; station < n_stations; ++station) { + for (size_t station = 0; station < n_stations; ++station) { // Make a list of all antennas for the current station, // and select the corresponding values from check_sum. std::vector<int> station_indices(n_receivers_per_station); std::vector<T> check_sum_station(n_receivers_per_station * 2); - for (int i = 0; i < n_receivers_per_station; ++i) { - int station_index = station * n_receivers_per_station + i; + for (size_t i = 0; i < n_receivers_per_station; ++i) { + const size_t station_index = station * n_receivers_per_station + i; station_indices[i] = station_index; T real = std::log(check_sum[station_index].real()); T imag = std::log(check_sum[station_index].imag()); @@ -175,7 +175,7 @@ void find_bad_antennas_(int n_stations, int n_receivers_per_station, } // Find outliers - int n_flags = n_receivers_per_station * 2; + const size_t n_flags = n_receivers_per_station * 2; std::vector<bool> flags = find_outliers(sigma, maxiters, check_sum_station); #if defined(DEBUG) @@ -184,7 +184,7 @@ void find_bad_antennas_(int n_stations, int n_receivers_per_station, // Count number of outliers int n_outliers = 0; - for (int i = 0; i < n_flags; ++i) { + for (size_t i = 0; i < n_flags; ++i) { n_outliers += flags[i]; } #if defined(DEBUG) @@ -196,12 +196,12 @@ void find_bad_antennas_(int n_stations, int n_receivers_per_station, #endif // Add the outliers to the list of bad stations - const int n = bad_stations[station].size(); + const size_t n = bad_stations[station].size(); bad_stations[station].resize(n + n_outliers); - int j = 0; - for (int i = 0; i < n_receivers_per_station * 2; ++i) { + size_t j = 0; + for (size_t i = 0; i < n_receivers_per_station * 2; ++i) { if (flags[i]) { - int station_index = + size_t station_index = station * n_receivers_per_station + (i % n_receivers_per_station); bad_stations[station][n + j++] = station_index; } @@ -210,7 +210,7 @@ void find_bad_antennas_(int n_stations, int n_receivers_per_station, } // Add the bad stations to the global list of bad stations - for (int station = 0; station < n_stations; ++station) { + for (size_t station = 0; station < n_stations; ++station) { all_bad_stations->insert(all_bad_stations->end(), bad_stations[station].begin(), bad_stations[station].end()); @@ -218,10 +218,10 @@ void find_bad_antennas_(int n_stations, int n_receivers_per_station, } // end find_bad_antennas_ template <typename T> -std::vector<std::complex<T>> compute_means(int n_stations, - int n_receivers_per_station, +std::vector<std::complex<T>> compute_means(size_t n_stations, + size_t n_receivers_per_station, const std::complex<T>* check_sum) { - const int n_receivers = n_stations * n_receivers_per_station; + const size_t n_receivers = n_stations * n_receivers_per_station; // - in: T2 check_sum[n_receivers] // - out: T2 means[n_stations] @@ -230,7 +230,7 @@ std::vector<std::complex<T>> compute_means(int n_stations, // Split the real and imaginary parts of the check sum std::vector<T> check_sum_real(n_receivers); std::vector<T> check_sum_imag(n_receivers); - for (int antenna = 0; antenna < n_receivers; ++antenna) { + for (size_t antenna = 0; antenna < n_receivers; ++antenna) { check_sum_real[antenna] = check_sum[antenna].real(); check_sum_imag[antenna] = check_sum[antenna].imag(); } @@ -239,13 +239,13 @@ std::vector<std::complex<T>> compute_means(int n_stations, T median_real = compute_median(check_sum_real); T median_imag = compute_median(check_sum_imag); - for (int station = 0; station < n_stations; ++station) { + for (size_t station = 0; station < n_stations; ++station) { // Make a list of all antennas for the current station, // and select the corresponding values from check_sum std::vector<T> check_sum_station_real(n_receivers_per_station); std::vector<T> check_sum_station_imag(n_receivers_per_station); - for (int i = 0; i < n_receivers_per_station; ++i) { - int station_index = station * n_receivers_per_station + i; + for (size_t i = 0; i < n_receivers_per_station; ++i) { + size_t station_index = station * n_receivers_per_station + i; check_sum_station_real[i] = check_sum_real[station_index]; check_sum_station_imag[i] = check_sum_imag[station_index]; } @@ -302,19 +302,19 @@ void AntennaFlagger::computeStats( } void AntennaFlagger::assertStatsComputed() const { - if (!(stats_std_.size() == size_t(n_correlations_ * n_receivers_)) || - !(stats_sump2_.size() == size_t(n_correlations_ * n_receivers_))) { + if (!(stats_std_.size() == (n_correlations_ * n_receivers_)) || + !(stats_sump2_.size() == (n_correlations_ * n_receivers_))) { throw std::runtime_error("Call AntennaFlagger::compute_stats() first!"); } } -std::vector<int> AntennaFlagger::findBadAntennas(float sigma, int maxiters, - int outlier_threshold) { +std::vector<size_t> AntennaFlagger::findBadAntennas(float sigma, int maxiters, + int outlier_threshold) { assertStatsComputed(); findbadAntennasWatch_.start(); - std::vector<int> flagged_antennas; + std::vector<size_t> flagged_antennas; find_bad_antennas_(n_stations_, n_receivers_per_station, n_baselines_, n_correlations_, sigma, maxiters, stats_sump2_, &flagged_antennas); @@ -323,15 +323,15 @@ std::vector<int> AntennaFlagger::findBadAntennas(float sigma, int maxiters, &flagged_antennas); // Count the number of times an antenna is flagged - std::vector<int> flagged_count(n_receivers_); - for (int antenna : flagged_antennas) { + std::vector<size_t> flagged_count(n_receivers_); + for (size_t antenna : flagged_antennas) { flagged_count[antenna]++; } // Make a list of antennas that are flagged multiple times - std::vector<int> flagged_antennas2; - for (int antenna = 0; antenna < n_receivers_; ++antenna) { - if (flagged_count[antenna] > outlier_threshold) { + std::vector<size_t> flagged_antennas2; + for (size_t antenna = 0; antenna < n_receivers_; ++antenna) { + if (flagged_count[antenna] > size_t(outlier_threshold)) { flagged_antennas2.push_back(antenna); } } @@ -341,7 +341,7 @@ std::vector<int> AntennaFlagger::findBadAntennas(float sigma, int maxiters, return flagged_antennas2; } -std::vector<int> AntennaFlagger::findBadStations(float sigma, int maxiters) { +std::vector<size_t> AntennaFlagger::findBadStations(float sigma, int maxiters) { assertStatsComputed(); findBadStationsWatch_.start(); @@ -362,7 +362,7 @@ std::vector<int> AntennaFlagger::findBadStations(float sigma, int maxiters) { for (size_t i = 0; i < check_sums.size(); ++i) { std::vector<std::complex<STATS_TYPE>> means = compute_means(n_stations_, n_receivers_per_station, check_sums[i]); - for (int station = 0; station < n_stations_; ++station) { + for (size_t station = 0; station < n_stations_; ++station) { all_means[0 * check_sums.size() + i][station] = means[station].real(); all_means[1 * check_sums.size() + i][station] = means[station].imag(); } @@ -379,23 +379,23 @@ std::vector<int> AntennaFlagger::findBadStations(float sigma, int maxiters) { assert(flags.size() == size_t(n_stations_)); #endif - for (int station = 0; station < n_stations_; ++station) { + for (size_t station = 0; station < n_stations_; ++station) { station_flags[station] = (flags[station] && (i == 0 || station_flags[station])); } } // Make a list of all the antenna indices for the flagged stations - std::vector<int> flagged_antennas; - for (int station = 0; station < n_stations_; ++station) { + std::vector<size_t> flagged_antennas; + for (size_t station = 0; station < n_stations_; ++station) { if (station_flags[station]) { // Make space for the new indices - int n = flagged_antennas.size(); + const size_t n = flagged_antennas.size(); flagged_antennas.resize(n + n_receivers_per_station); // Add the new indices - for (int i = 0; i < n_receivers_per_station; ++i) { - int station_index = station * n_receivers_per_station + i; + for (size_t i = 0; i < n_receivers_per_station; ++i) { + const size_t station_index = station * n_receivers_per_station + i; flagged_antennas[n + i] = station_index; } } diff --git a/antennaflagger/AntennaFlagger.h b/antennaflagger/AntennaFlagger.h index 26dec964d..484ce9bd9 100644 --- a/antennaflagger/AntennaFlagger.h +++ b/antennaflagger/AntennaFlagger.h @@ -28,8 +28,8 @@ namespace antennaflagger { */ class AntennaFlagger { public: - AntennaFlagger(int n_stations, int n_receivers_per_station, int n_channels, - int n_correlations) + AntennaFlagger(size_t n_stations, size_t n_receivers_per_station, + size_t n_channels, size_t n_correlations) : n_stations_(n_stations), n_receivers_per_station(n_receivers_per_station), n_channels_(n_channels), @@ -40,9 +40,9 @@ class AntennaFlagger { void computeStats(const std::vector<std::complex<double>>& data); void computeStats(const std::vector<std::complex<float>>& data); - std::vector<int> findBadAntennas(float sigma, int maxiters, - int outlier_threshold); - std::vector<int> findBadStations(float sigma, int maxiters); + std::vector<size_t> findBadAntennas(float sigma, int maxiters, + int outlier_threshold); + std::vector<size_t> findBadStations(float sigma, int maxiters); void reportRuntime() const; @@ -51,14 +51,14 @@ class AntennaFlagger { void assertStatsComputed() const; // Constants - const int n_stations_; - const int n_receivers_per_station; - const int n_channels_; - const int n_correlations_; + const size_t n_stations_; + const size_t n_receivers_per_station; + const size_t n_channels_; + const size_t n_correlations_; // Derived constants - const int n_receivers_; - const int n_baselines_; + const size_t n_receivers_; + const size_t n_baselines_; // Timings common::NSTimer computeStatisticsWatch_; diff --git a/common/baseline_indices.h b/common/baseline_indices.h index 7f0a9b13f..763d6e232 100644 --- a/common/baseline_indices.h +++ b/common/baseline_indices.h @@ -65,11 +65,10 @@ enum class BaselineOrder { kColumnMajor, kRowMajor }; /* * Computes the number of baselines, including auto-correlations. */ -inline int ComputeNBaselines(int n_antennas) { - int n_correlations = (n_antennas * (n_antennas - 1)) / 2; - int n_autocorrelations = n_antennas; - int n_baselines = n_correlations + n_autocorrelations; - return n_baselines; +inline size_t ComputeNBaselines(size_t n_antennas) { + const size_t n_correlations = (n_antennas * (n_antennas - 1)) / 2; + const size_t n_autocorrelations = n_antennas; + return n_correlations + n_autocorrelations; } /* @@ -77,14 +76,14 @@ inline int ComputeNBaselines(int n_antennas) { * number of antennas. It assumes that auto-correlations of antennas are * included. */ -inline int ComputeBaselineIndex(int antenna_a, int antenna_b, int nAntennas, - BaselineOrder order) { +inline size_t ComputeBaselineIndex(size_t antenna_a, size_t antenna_b, + size_t nAntennas, BaselineOrder order) { size_t row = std::min(antenna_a, antenna_b); size_t col = std::max(antenna_a, antenna_b); if (order == BaselineOrder::kRowMajor) { - return int((row * nAntennas) + col - row * (row + 1) / 2); + return size_t((row * nAntennas) + col - row * (row + 1) / 2); } else { - return int((col * (col + 1)) / 2) + row; + return size_t((col * (col + 1)) / 2) + row; } } @@ -92,22 +91,23 @@ inline int ComputeBaselineIndex(int antenna_a, int antenna_b, int nAntennas, * Computes the antenna pair (baseline) given the baseline index. It assumes * that auto-correlations of antennas are included. */ -inline std::pair<int, int> ComputeBaseline(int baseline_index, int n_antennas, - BaselineOrder order) { +inline std::pair<size_t, size_t> ComputeBaseline(size_t baseline_index, + size_t n_antennas, + BaselineOrder order) { if (order == BaselineOrder::kRowMajor) { - const int n_baselines = ComputeNBaselines(n_antennas); - const int antenna = + const size_t n_baselines = ComputeNBaselines(n_antennas); + const size_t antenna = (1 + std::sqrt(1 + 8 * (n_baselines - baseline_index - 1))) / 2; - const int antenna1 = n_antennas - antenna; + const size_t antenna1 = n_antennas - antenna; // n is the number of baselines (a,b) with a < antenna - const int n = n_baselines - n_antennas - (antenna * (antenna - 1)) / 2; - const int antenna2 = baseline_index - n; + const size_t n = n_baselines - n_antennas - (antenna * (antenna - 1)) / 2; + const size_t antenna2 = baseline_index - n; return {antenna1, antenna2}; } else { - const int antenna2 = (1 + std::sqrt(1 + 8 * baseline_index)) / 2; + const size_t antenna2 = (1 + std::sqrt(1 + 8 * baseline_index)) / 2; // n is the number of baselines (a,b) with b < antenna2 - const int n = int((antenna2 * (antenna2 - 1)) / 2); - const int antenna1 = baseline_index - n; + const size_t n = int((antenna2 * (antenna2 - 1)) / 2); + const size_t antenna1 = baseline_index - n; return {antenna1, antenna2 - 1}; } } @@ -116,12 +116,13 @@ inline std::pair<int, int> ComputeBaseline(int baseline_index, int n_antennas, * Computes the list of baselines for a given antenna. It assumes that * auto-correlations are included. */ -inline std::vector<int> ComputeBaselineList(int antenna, int n_antennas, - BaselineOrder order) { - std::vector<int> indices; +inline std::vector<size_t> ComputeBaselineList(size_t antenna, + size_t n_antennas, + BaselineOrder order) { + std::vector<size_t> indices; indices.reserve(n_antennas); - for (int i = 0; i < n_antennas; i++) { + for (size_t i = 0; i < n_antennas; i++) { indices.push_back(ComputeBaselineIndex(antenna, i, n_antennas, order)); } diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index dce1892ab..d72642dbe 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -207,15 +207,15 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { flagger_->computeStats(data); // Find outliers both at antenna and station level - const std::vector<int> flagged_antennas_1 = flagger_->findBadAntennas( + const std::vector<size_t> flagged_antennas_1 = flagger_->findBadAntennas( antenna_flagging_sigma_, antenna_flagging_maxiters_, antenna_flagging_outlier_threshold_); - const std::vector<int> flagged_antennas_2 = flagger_->findBadStations( + const std::vector<size_t> flagged_antennas_2 = flagger_->findBadStations( station_flagging_sigma_, station_flagging_maxiters_); // Combine vectors into a single set. - std::set<int> flagged_antennas(flagged_antennas_2.begin(), - flagged_antennas_2.end()); + std::set<size_t> flagged_antennas(flagged_antennas_2.begin(), + flagged_antennas_2.end()); flagged_antennas.insert(flagged_antennas_1.begin(), flagged_antennas_1.end()); -- GitLab From 73f1744151f1d27a47d94d95640f7eb2fea43a9d Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 08:50:57 +0200 Subject: [PATCH 39/61] Update AntennaFlagger constructor --- steps/AntennaFlagger.cc | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index d72642dbe..81e7f27a7 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -162,23 +162,21 @@ AntennaFlagger::AntennaFlagger(InputStep* input, selection_ = Convert(selectionString, n_antennas_); name_ = prefix; - // Initialize flagger class if (doDetectOutliers_) { - // Find out if the observation has LBA, HBA or AARTFAAC-12 data. std::string antennaSet(input->getInfo().antennaSet()); std::vector<std::string> antennaNames(input->getInfo().antennaNames()); size_t n_receivers_per_station = 0; size_t n_stations = 0; - if (antennaSet.substr(0, 3) == "LBA" || antennaSet.substr(0, 3) == "HBA") { - // LOFAR LBA or HBA - n_receivers_per_station = 1; - n_stations = input->getInfo().nantenna(); - } else if (antennaNames[0].substr(0, 3) == "A12") { - // AARTFAAC-12 + if (antennaNames[0].substr(0, 3) == "A12") { + // AARTFAAC-12: all antennas in a station are used as individual receivers n_receivers_per_station = 48; - n_stations = 12; + const size_t n_antennas = input->getInfo().nantenna(); + n_stations = n_antennas / n_receivers_per_station; } else { - throw std::runtime_error("Could not determine observation type."); + // e.g. LOFAR LBA or HBA: all antennas in a station are combined to form + // one 'receiver' + n_receivers_per_station = 1; + n_stations = input->getInfo().nantenna(); } flagger_ = std::make_unique<dp3::antennaflagger::AntennaFlagger>( -- GitLab From f24ccbe3cdb0917b0ed4f32dbe17411ff9e3e6a6 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 08:55:07 +0200 Subject: [PATCH 40/61] Cleanup debug code --- antennaflagger/statistics.h | 17 +---------------- 1 file changed, 1 insertion(+), 16 deletions(-) diff --git a/antennaflagger/statistics.h b/antennaflagger/statistics.h index be7a25d55..909353315 100644 --- a/antennaflagger/statistics.h +++ b/antennaflagger/statistics.h @@ -1,18 +1,12 @@ #include <algorithm> -#include <vector> -#if defined(DEBUG) #include <cassert> -#include <iostream> -#include <iomanip> -#endif +#include <vector> namespace { template <typename T> T compute_mean(const std::vector<T>& data, const std::vector<bool>& flags) { -#if defined(DEBUG) assert(data.size() == flags.size()); -#endif T sum = 0; size_t m = 0; @@ -124,9 +118,7 @@ T compute_median(std::vector<T> v) { template <typename T> T compute_std(const std::vector<T>& data, const std::vector<bool>& flags) { -#if defined(DEBUG) assert(data.size() == flags.size()); -#endif T sum = 0; size_t m = 0; @@ -174,9 +166,7 @@ inline bool is_valid(const std::complex<T>& x) { template <typename T> void flag_invalid(const std::vector<T>& data, std::vector<bool>& flags) { -#if defined(DEBUG) assert(data.size() == flags.size()); -#endif for (size_t i = 0; i < data.size(); i++) { if (!is_valid(data[i])) { @@ -208,11 +198,6 @@ std::vector<bool> find_outliers(float sigma, int maxiters, if (!flags[j] && (value < lower_bound || value > upper_bound)) { flags[j] = true; m++; -#if defined(DEBUG) - std::cout << "flagging data[" << std::setw(2) << std::scientific << j - << "] = " << value << " (bounds: [" << lower_bound << "," - << upper_bound << "]" << std::endl; -#endif } } -- GitLab From 83d0aac130b8e9cde8caa052082d3c3a29763ef8 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 10:31:14 +0200 Subject: [PATCH 41/61] Update statistics --- antennaflagger/AntennaFlagger.cc | 12 +++--- antennaflagger/statistics.h | 72 ++++++++++++++++---------------- 2 files changed, 42 insertions(+), 42 deletions(-) diff --git a/antennaflagger/AntennaFlagger.cc b/antennaflagger/AntennaFlagger.cc index 16a42c4f1..e8a5aaab2 100644 --- a/antennaflagger/AntennaFlagger.cc +++ b/antennaflagger/AntennaFlagger.cc @@ -177,7 +177,7 @@ void find_bad_antennas_(size_t n_stations, size_t n_receivers_per_station, // Find outliers const size_t n_flags = n_receivers_per_station * 2; std::vector<bool> flags = - find_outliers(sigma, maxiters, check_sum_station); + FindOutliers(sigma, maxiters, check_sum_station); #if defined(DEBUG) assert(flags.size() == size_t(n_flags)); #endif @@ -236,8 +236,8 @@ std::vector<std::complex<T>> compute_means(size_t n_stations, } // Compute the median for both real and imaginary parts - T median_real = compute_median(check_sum_real); - T median_imag = compute_median(check_sum_imag); + T median_real = ComputeMedian(check_sum_real); + T median_imag = ComputeMedian(check_sum_imag); for (size_t station = 0; station < n_stations; ++station) { // Make a list of all antennas for the current station, @@ -252,8 +252,8 @@ std::vector<std::complex<T>> compute_means(size_t n_stations, // Compute the median for both the real and imaginary parts // of the check sum for the current station - T median_station_real = compute_median(check_sum_station_real); - T median_station_imag = compute_median(check_sum_station_imag); + T median_station_real = ComputeMedian(check_sum_station_real); + T median_station_imag = ComputeMedian(check_sum_station_imag); // Add the mean value for the current station means[station] = {median_station_real / median_real, @@ -374,7 +374,7 @@ std::vector<size_t> AntennaFlagger::findBadStations(float sigma, int maxiters) { const std::vector<STATS_TYPE>& means = all_means[i]; // Find outliers - std::vector<bool> flags = find_outliers(sigma, maxiters, means); + std::vector<bool> flags = FindOutliers(sigma, maxiters, means); #if defined(DEBUG) assert(flags.size() == size_t(n_stations_)); #endif diff --git a/antennaflagger/statistics.h b/antennaflagger/statistics.h index 909353315..41c851505 100644 --- a/antennaflagger/statistics.h +++ b/antennaflagger/statistics.h @@ -5,13 +5,13 @@ namespace { template <typename T> -T compute_mean(const std::vector<T>& data, const std::vector<bool>& flags) { +T ComputeMean(const std::vector<T>& data, const std::vector<bool>& flags) { assert(data.size() == flags.size()); T sum = 0; size_t m = 0; - for (size_t i = 0; i < data.size(); i++) { + for (size_t i = 0; i < data.size(); ++i) { if (!flags[i]) { sum += data[i]; m++; @@ -29,28 +29,28 @@ T compute_mean(const std::vector<T>& data, const std::vector<bool>& flags) { * with a number of notable changes: * - The input is an array with values of type T (instead of objects with value * and weight) - * - Skips flagged values (NaN or Inf) - * - Handles even-length input as a special case + * - Skips values for which its flag is true, or which are NaN or Inf + * - Handles even-length input as a special case */ template <typename T> -T compute_median_inplace(const std::vector<T>& data, - const std::vector<bool>& flags) { +T ComputeMedianInplace(const std::vector<T>& data, + const std::vector<bool>& flags) { assert(data.size() == flags.size()); T below = 0; T above = 0; - T equal = compute_mean(data, flags); + T equal = ComputeMean(data, flags); while (true) { size_t n_below = 0; size_t n_equal = 0; size_t n_above = 0; - for (size_t i = 0; i < data.size(); i++) { - if (flags[i]) { + for (size_t i = 0; i < data.size(); ++i) { + T value = data[i]; + if (flags[i] || !IsValid(value)) { continue; } - T value = data[i]; if (value > equal) { if (n_above == 0 || value < above) { above = value; @@ -69,7 +69,7 @@ T compute_median_inplace(const std::vector<T>& data, T mid = n_below + n_equal + n_above; if (mid > 0) { - mid /= 2; + mid /= static_cast<T>(2); } else { return 0; } @@ -80,7 +80,7 @@ T compute_median_inplace(const std::vector<T>& data, return equal; } else { if (mid == n_below) { - for (size_t j = 0; j < data.size(); j++) { + for (size_t j = 0; j < data.size(); ++j) { if (flags[j]) { continue; } @@ -89,9 +89,9 @@ T compute_median_inplace(const std::vector<T>& data, above = value; } } - return (below + above) / 2; + return (below + above) / static_cast<T>(2); } else { - return (equal + above) / 2; + return (equal + above) / static_cast<T>(2); } } } @@ -104,8 +104,8 @@ T compute_median_inplace(const std::vector<T>& data, } template <typename T> -T compute_median(std::vector<T> v) { - auto target = v.begin() + v.size() / 2; +T ComputeMedian(std::vector<T> v) { + auto target = v.begin() + v.size() / static_cast<T>(2); std::nth_element(v.begin(), target, v.end()); if (v.size() % 2 != 0) { return *target; @@ -117,15 +117,15 @@ T compute_median(std::vector<T> v) { } template <typename T> -T compute_std(const std::vector<T>& data, const std::vector<bool>& flags) { +T ComputeSTD(const std::vector<T>& data, const std::vector<bool>& flags) { assert(data.size() == flags.size()); T sum = 0; size_t m = 0; - T mean = compute_mean(data, flags); + T mean = ComputeMean(data, flags); - for (size_t i = 0; i < data.size(); i++) { + for (size_t i = 0; i < data.size(); ++i) { if (!flags[i]) { T value = mean - data[i]; sum += value * value; @@ -133,16 +133,16 @@ T compute_std(const std::vector<T>& data, const std::vector<bool>& flags) { } } - return m > 0 ? std::sqrt(sum / m) : 0; + return m > 0 ? std::sqrt(sum / m) : static_cast<T>(0); } template <typename T> -T compute_sump2(const std::vector<T>& data, const std::vector<bool>& flags) { +T ComputeSUMP2(const std::vector<T>& data, const std::vector<bool>& flags) { assert(data.size() == flags.size()); T sum = 0; - for (size_t i = 0; i < data.size(); i++) { + for (size_t i = 0; i < data.size(); ++i) { if (!flags[i]) { T value = data[i]; sum += value * value; @@ -153,48 +153,48 @@ T compute_sump2(const std::vector<T>& data, const std::vector<bool>& flags) { } template <typename T> -inline bool is_valid(T x) { +inline bool IsValid(T x) { bool is_nan = x != x; bool is_finite = x != 2 * x; return !is_nan && is_finite; } template <typename T> -inline bool is_valid(const std::complex<T>& x) { - return is_valid(x.real()) && is_valid(x.imag()); +inline bool IsValid(const std::complex<T>& x) { + return IsValid(x.real()) && IsValid(x.imag()); } template <typename T> -void flag_invalid(const std::vector<T>& data, std::vector<bool>& flags) { +void FlagInvalid(const std::vector<T>& data, std::vector<bool>& flags) { assert(data.size() == flags.size()); - for (size_t i = 0; i < data.size(); i++) { - if (!is_valid(data[i])) { + for (size_t i = 0; i < data.size(); ++i) { + if (!IsValid(data[i])) { flags[i] = true; } } } template <typename T> -std::vector<bool> find_outliers(float sigma, int maxiters, - const std::vector<T>& data) { +std::vector<bool> FindOutliers(float sigma, int maxiters, + const std::vector<T>& data) { std::vector<bool> flags(data.size(), {false}); - flag_invalid(data, flags); + FlagInvalid(data, flags); - for (int i = 0; i < maxiters; i++) { + for (int i = 0; i < maxiters; ++i) { // Compute median - T median = compute_median(data); + T median = ComputeMedian(data); // Compute standard deviation - T std = compute_std(data, flags); + T std = ComputeSTD(data, flags); // Find outliers size_t m = 0; T lower_bound = median - (sigma * std); T upper_bound = median + (sigma * std); - for (size_t j = 0; j < data.size(); j++) { - T value = data[j]; + for (size_t j = 0; j < data.size(); ++j) { + const T value = data[j]; if (!flags[j] && (value < lower_bound || value > upper_bound)) { flags[j] = true; m++; -- GitLab From 784b7ecff50e4ad092fb79ea330868210b44557c Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 10:33:00 +0200 Subject: [PATCH 42/61] Remove anonymous namespace, add namespace detail --- antennaflagger/statistics.h | 56 ++++++++++++++++++------------------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/antennaflagger/statistics.h b/antennaflagger/statistics.h index 41c851505..60e6a5d45 100644 --- a/antennaflagger/statistics.h +++ b/antennaflagger/statistics.h @@ -2,7 +2,32 @@ #include <cassert> #include <vector> -namespace { +namespace detail { + +template <typename T> +inline bool IsValid(T x) { + bool is_nan = x != x; + bool is_finite = x != 2 * x; + return !is_nan && is_finite; +} + +template <typename T> +inline bool IsValid(const std::complex<T>& x) { + return IsValid(x.real()) && IsValid(x.imag()); +} + +template <typename T> +void FlagInvalid(const std::vector<T>& data, std::vector<bool>& flags) { + assert(data.size() == flags.size()); + + for (size_t i = 0; i < data.size(); ++i) { + if (!IsValid(data[i])) { + flags[i] = true; + } + } +} + +} // namespace detail template <typename T> T ComputeMean(const std::vector<T>& data, const std::vector<bool>& flags) { @@ -152,34 +177,11 @@ T ComputeSUMP2(const std::vector<T>& data, const std::vector<bool>& flags) { return sum; } -template <typename T> -inline bool IsValid(T x) { - bool is_nan = x != x; - bool is_finite = x != 2 * x; - return !is_nan && is_finite; -} - -template <typename T> -inline bool IsValid(const std::complex<T>& x) { - return IsValid(x.real()) && IsValid(x.imag()); -} - -template <typename T> -void FlagInvalid(const std::vector<T>& data, std::vector<bool>& flags) { - assert(data.size() == flags.size()); - - for (size_t i = 0; i < data.size(); ++i) { - if (!IsValid(data[i])) { - flags[i] = true; - } - } -} - template <typename T> std::vector<bool> FindOutliers(float sigma, int maxiters, const std::vector<T>& data) { std::vector<bool> flags(data.size(), {false}); - FlagInvalid(data, flags); + detail::FlagInvalid(data, flags); for (int i = 0; i < maxiters; ++i) { // Compute median @@ -208,6 +210,4 @@ std::vector<bool> FindOutliers(float sigma, int maxiters, } return flags; -} - -} // anonymous namespace \ No newline at end of file +} \ No newline at end of file -- GitLab From d87f54d3fd10fbda7ad39b7c006ea5efab0122e2 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 10:34:41 +0200 Subject: [PATCH 43/61] Add missing include --- antennaflagger/statistics.h | 1 + 1 file changed, 1 insertion(+) diff --git a/antennaflagger/statistics.h b/antennaflagger/statistics.h index 60e6a5d45..3f5f74aba 100644 --- a/antennaflagger/statistics.h +++ b/antennaflagger/statistics.h @@ -1,5 +1,6 @@ #include <algorithm> #include <cassert> +#include <complex> #include <vector> namespace detail { -- GitLab From 8ecea384d6509276c7d352ee5e75cd5eff2e6397 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 10:39:29 +0200 Subject: [PATCH 44/61] Update statistics --- antennaflagger/statistics.h | 33 +++++++++++++++------------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/antennaflagger/statistics.h b/antennaflagger/statistics.h index 3f5f74aba..f95376208 100644 --- a/antennaflagger/statistics.h +++ b/antennaflagger/statistics.h @@ -35,16 +35,16 @@ T ComputeMean(const std::vector<T>& data, const std::vector<bool>& flags) { assert(data.size() == flags.size()); T sum = 0; - size_t m = 0; + size_t unflagged_count = 0; for (size_t i = 0; i < data.size(); ++i) { if (!flags[i]) { sum += data[i]; - m++; + ++unflagged_count; } } - return m > 0 ? sum / m : 0; + return unflagged_count > 0 ? sum / unflagged_count : 0; } /* @@ -81,14 +81,14 @@ T ComputeMedianInplace(const std::vector<T>& data, if (n_above == 0 || value < above) { above = value; } - n_above++; + ++n_above; } else if (value < equal) { if (n_below == 0 || value > below) { below = value; } - n_below++; + ++n_below; } else { - n_equal++; + ++n_equal; } } @@ -155,7 +155,7 @@ T ComputeSTD(const std::vector<T>& data, const std::vector<bool>& flags) { if (!flags[i]) { T value = mean - data[i]; sum += value * value; - m++; + ++m; } } @@ -170,7 +170,7 @@ T ComputeSUMP2(const std::vector<T>& data, const std::vector<bool>& flags) { for (size_t i = 0; i < data.size(); ++i) { if (!flags[i]) { - T value = data[i]; + const T value = data[i]; sum += value * value; } } @@ -185,27 +185,24 @@ std::vector<bool> FindOutliers(float sigma, int maxiters, detail::FlagInvalid(data, flags); for (int i = 0; i < maxiters; ++i) { - // Compute median - T median = ComputeMedian(data); - - // Compute standard deviation - T std = ComputeSTD(data, flags); + const T median = ComputeMedian(data); + const T std = ComputeSTD(data, flags); // Find outliers - size_t m = 0; - T lower_bound = median - (sigma * std); - T upper_bound = median + (sigma * std); + size_t outlier_count = 0; + const T lower_bound = median - (sigma * std); + const T upper_bound = median + (sigma * std); for (size_t j = 0; j < data.size(); ++j) { const T value = data[j]; if (!flags[j] && (value < lower_bound || value > upper_bound)) { flags[j] = true; - m++; + ++outlier_count; } } // Exit early when no new outliers were found - if (m == 0) { + if (outlier_count == 0) { break; } } -- GitLab From 761ffd27e84e29988363b6e259a6477f5fd3f635 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 10:45:54 +0200 Subject: [PATCH 45/61] Rename AntennaFlagger implementation to Flagger --- CMakeLists.txt | 2 +- .../{AntennaFlagger.cc => Flagger.cc} | 20 +++++++++---------- .../{AntennaFlagger.h => Flagger.h} | 6 +++--- steps/AntennaFlagger.cc | 2 +- steps/AntennaFlagger.h | 4 ++-- 5 files changed, 16 insertions(+), 18 deletions(-) rename antennaflagger/{AntennaFlagger.cc => Flagger.cc} (96%) rename antennaflagger/{AntennaFlagger.h => Flagger.h} (94%) diff --git a/CMakeLists.txt b/CMakeLists.txt index b3f992604..8c20bcdb5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -382,7 +382,7 @@ add_library( common/VdsMaker.cc common/VdsPartDesc.cc) -add_library(AntennaFlagger OBJECT antennaflagger/AntennaFlagger.cc) +add_library(AntennaFlagger OBJECT antennaflagger/Flagger.cc) set_property(TARGET AntennaFlagger PROPERTY POSITION_INDEPENDENT_CODE ON) add_library( diff --git a/antennaflagger/AntennaFlagger.cc b/antennaflagger/Flagger.cc similarity index 96% rename from antennaflagger/AntennaFlagger.cc rename to antennaflagger/Flagger.cc index e8a5aaab2..b4b62a5a5 100644 --- a/antennaflagger/AntennaFlagger.cc +++ b/antennaflagger/Flagger.cc @@ -8,7 +8,7 @@ #include <iostream> #include <set> -#include "AntennaFlagger.h" +#include "Flagger.h" #include "./statistics.h" // Non-member helper functions @@ -271,14 +271,13 @@ std::vector<std::complex<T>> compute_means(size_t n_stations, namespace dp3 { namespace antennaflagger { -void AntennaFlagger::reportRuntime() const { +void Flagger::reportRuntime() const { std::cout << "statistics: " << computeStatisticsWatch_.getElapsed() << ", bad_antennas: " << findbadAntennasWatch_.getElapsed() << ", bad_stations: " << findBadStationsWatch_.getElapsed() << '\n'; } -void AntennaFlagger::computeStats( - const std::vector<std::complex<double>>& data) { +void Flagger::computeStats(const std::vector<std::complex<double>>& data) { computeStatisticsWatch_.start(); compute_stats_per_antenna(n_receivers_, n_baselines_, n_channels_, n_correlations_, data, &stats_std_, &stats_sump2_); @@ -289,8 +288,7 @@ void AntennaFlagger::computeStats( computeStatisticsWatch_.stop(); } -void AntennaFlagger::computeStats( - const std::vector<std::complex<float>>& data) { +void Flagger::computeStats(const std::vector<std::complex<float>>& data) { computeStatisticsWatch_.start(); compute_stats_per_antenna(n_receivers_, n_baselines_, n_channels_, n_correlations_, data, &stats_std_, &stats_sump2_); @@ -301,15 +299,15 @@ void AntennaFlagger::computeStats( computeStatisticsWatch_.stop(); } -void AntennaFlagger::assertStatsComputed() const { +void Flagger::assertStatsComputed() const { if (!(stats_std_.size() == (n_correlations_ * n_receivers_)) || !(stats_sump2_.size() == (n_correlations_ * n_receivers_))) { - throw std::runtime_error("Call AntennaFlagger::compute_stats() first!"); + throw std::runtime_error("Call Flagger::compute_stats() first!"); } } -std::vector<size_t> AntennaFlagger::findBadAntennas(float sigma, int maxiters, - int outlier_threshold) { +std::vector<size_t> Flagger::findBadAntennas(float sigma, int maxiters, + int outlier_threshold) { assertStatsComputed(); findbadAntennasWatch_.start(); @@ -341,7 +339,7 @@ std::vector<size_t> AntennaFlagger::findBadAntennas(float sigma, int maxiters, return flagged_antennas2; } -std::vector<size_t> AntennaFlagger::findBadStations(float sigma, int maxiters) { +std::vector<size_t> Flagger::findBadStations(float sigma, int maxiters) { assertStatsComputed(); findBadStationsWatch_.start(); diff --git a/antennaflagger/AntennaFlagger.h b/antennaflagger/Flagger.h similarity index 94% rename from antennaflagger/AntennaFlagger.h rename to antennaflagger/Flagger.h index 484ce9bd9..1682a0f80 100644 --- a/antennaflagger/AntennaFlagger.h +++ b/antennaflagger/Flagger.h @@ -26,10 +26,10 @@ namespace antennaflagger { * 3) Query these statistics to find outlier antennas * 4) Query these statistics to find outlier stations */ -class AntennaFlagger { +class Flagger { public: - AntennaFlagger(size_t n_stations, size_t n_receivers_per_station, - size_t n_channels, size_t n_correlations) + Flagger(size_t n_stations, size_t n_receivers_per_station, size_t n_channels, + size_t n_correlations) : n_stations_(n_stations), n_receivers_per_station(n_receivers_per_station), n_channels_(n_channels), diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index 81e7f27a7..c6cee7af6 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -179,7 +179,7 @@ AntennaFlagger::AntennaFlagger(InputStep* input, n_stations = input->getInfo().nantenna(); } - flagger_ = std::make_unique<dp3::antennaflagger::AntennaFlagger>( + flagger_ = std::make_unique<dp3::antennaflagger::Flagger>( n_stations, n_receivers_per_station, n_channels_, n_correlations_); } diff --git a/steps/AntennaFlagger.h b/steps/AntennaFlagger.h index d309be24a..7bd5048bd 100644 --- a/steps/AntennaFlagger.h +++ b/steps/AntennaFlagger.h @@ -9,7 +9,7 @@ #include <vector> #include <xtensor.hpp> -#include "../antennaflagger/AntennaFlagger.h" +#include "../antennaflagger/Flagger.h" #include "InputStep.h" @@ -37,7 +37,7 @@ class AntennaFlagger final : public Step { std::string name_; std::string selectionString_; bool doDetectOutliers_; - std::unique_ptr<dp3::antennaflagger::AntennaFlagger> flagger_; + std::unique_ptr<dp3::antennaflagger::Flagger> flagger_; // AartfaacFlagger parameters float antenna_flagging_sigma_; -- GitLab From 2a0a54ae5ad3b6cd18b6c411572d1e3bd46de533 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 10:54:03 +0200 Subject: [PATCH 46/61] Update Flagger --- antennaflagger/Flagger.cc | 93 ++++++++++++++++++++------------------- antennaflagger/Flagger.h | 12 ++--- steps/AntennaFlagger.cc | 6 +-- 3 files changed, 56 insertions(+), 55 deletions(-) diff --git a/antennaflagger/Flagger.cc b/antennaflagger/Flagger.cc index b4b62a5a5..e37c39320 100644 --- a/antennaflagger/Flagger.cc +++ b/antennaflagger/Flagger.cc @@ -11,10 +11,9 @@ #include "Flagger.h" #include "./statistics.h" -// Non-member helper functions -namespace { -inline size_t data_index(size_t n_channels, size_t n_correlations, size_t bl, - size_t chan, size_t cor) { +namespace detail { +inline size_t DataIndex(size_t n_channels, size_t n_correlations, size_t bl, + size_t chan, size_t cor) { // data[n_baselines][n_channels][n_correlations] return (bl * n_channels + chan) * n_correlations + cor; } @@ -23,11 +22,11 @@ inline size_t data_index(size_t n_channels, size_t n_correlations, size_t bl, * Helper functions for AartfaacFlagger::find_bad_antennas */ template <typename InType, typename OutType> -void compute_stats_per_baseline( - size_t n_baselines, size_t n_channels, size_t n_correlations, - const std::vector<std::complex<InType>>& data, - std::vector<std::complex<OutType>>* stats_std, - std::vector<std::complex<OutType>>* stats_sump2) { +void ComputeStatsPerBaseline(size_t n_baselines, size_t n_channels, + size_t n_correlations, + const std::vector<std::complex<InType>>& data, + std::vector<std::complex<OutType>>* stats_std, + std::vector<std::complex<OutType>>* stats_sump2) { // - in: data[n_baselines][n_channels][n_correlations] // - out: stats[n_correlations][n_baselines] // Note that the baseline and correlation indices are swapped between @@ -52,7 +51,7 @@ void compute_stats_per_baseline( OutType sump2_real = 0; OutType sump2_imag = 0; for (size_t chan = 0; chan < n_channels; ++chan) { - size_t idx = data_index(n_channels, n_correlations, bl, chan, cor); + size_t idx = DataIndex(n_channels, n_correlations, bl, chan, cor); std::complex<InType> value = data[idx]; sum_real += value.real(); sum_imag += value.imag(); @@ -68,7 +67,7 @@ void compute_stats_per_baseline( OutType sum2_real = 0; OutType sum2_imag = 0; for (size_t chan = 0; chan < n_channels; ++chan) { - size_t idx = data_index(n_channels, n_correlations, bl, chan, cor); + size_t idx = DataIndex(n_channels, n_correlations, bl, chan, cor); std::complex<InType> value = data[idx]; OutType diff_real = mean_real - value.real(); OutType diff_imag = mean_imag - value.imag(); @@ -84,10 +83,10 @@ void compute_stats_per_baseline( (*stats_sump2)[idx] = std::complex<OutType>(sump2_real, sump2_imag); } } -} // end compute_stats_per_baseline +} // end ComputeStatsPerBaseline template <typename InType, typename OutType> -void compute_stats_per_antenna( +void ComputeStatsPerAntenna( size_t n_receivers, size_t n_baselines, size_t n_channels, size_t n_correlations, const std::vector<std::complex<InType>>& data, std::vector<std::complex<OutType>>* stats_std_grouped, @@ -101,8 +100,8 @@ void compute_stats_per_antenna( #endif std::vector<std::complex<OutType>> stats_std; std::vector<std::complex<OutType>> stats_sump2; - compute_stats_per_baseline(n_baselines, n_channels, n_correlations, data, - &stats_std, &stats_sump2); + ComputeStatsPerBaseline(n_baselines, n_channels, n_correlations, data, + &stats_std, &stats_sump2); #if defined(DEBUG) assert(stats_std.size() == size_t(n_correlations * n_baselines)); assert(stats_sump2.size() == size_t(n_correlations * n_baselines)); @@ -131,16 +130,16 @@ void compute_stats_per_antenna( } } } -} // end compute_stats_per_antenna +} // end ComputeStatsPerAntenna /* * Helper functions for AaatfaacFlagger::find_bad_stations */ template <typename T> -void find_bad_antennas_(size_t n_stations, size_t n_receivers_per_station, - size_t n_baselines, size_t n_correlations, float sigma, - int maxiters, const std::vector<std::complex<T>>& stats, - std::vector<size_t>* all_bad_stations) { +void FindBadAntennas(size_t n_stations, size_t n_receivers_per_station, + size_t n_baselines, size_t n_correlations, float sigma, + int maxiters, const std::vector<std::complex<T>>& stats, + std::vector<size_t>* all_bad_stations) { int n_receivers = n_stations * n_receivers_per_station; // - in: stats[n_correlations][n_receivers] // - out: all_bad_stations[?] @@ -215,12 +214,12 @@ void find_bad_antennas_(size_t n_stations, size_t n_receivers_per_station, bad_stations[station].begin(), bad_stations[station].end()); } -} // end find_bad_antennas_ +} // end FindBadAntennas template <typename T> -std::vector<std::complex<T>> compute_means(size_t n_stations, - size_t n_receivers_per_station, - const std::complex<T>* check_sum) { +std::vector<std::complex<T>> ComputeMeans(size_t n_stations, + size_t n_receivers_per_station, + const std::complex<T>* check_sum) { const size_t n_receivers = n_stations * n_receivers_per_station; // - in: T2 check_sum[n_receivers] @@ -261,9 +260,9 @@ std::vector<std::complex<T>> compute_means(size_t n_stations, } return means; -} // end compute_means +} // end ComputeMeans -}; // namespace +}; // namespace detail /* * Member functions @@ -271,16 +270,17 @@ std::vector<std::complex<T>> compute_means(size_t n_stations, namespace dp3 { namespace antennaflagger { -void Flagger::reportRuntime() const { +void Flagger::ReportRuntime() const { std::cout << "statistics: " << computeStatisticsWatch_.getElapsed() << ", bad_antennas: " << findbadAntennasWatch_.getElapsed() << ", bad_stations: " << findBadStationsWatch_.getElapsed() << '\n'; } -void Flagger::computeStats(const std::vector<std::complex<double>>& data) { +void Flagger::ComputeStats(const std::vector<std::complex<double>>& data) { computeStatisticsWatch_.start(); - compute_stats_per_antenna(n_receivers_, n_baselines_, n_channels_, - n_correlations_, data, &stats_std_, &stats_sump2_); + detail::ComputeStatsPerAntenna(n_receivers_, n_baselines_, n_channels_, + n_correlations_, data, &stats_std_, + &stats_sump2_); #if defined(DEBUG) assert(stats_std_.size() == size_t(n_correlations_ * n_receivers_)); assert(stats_sump2_.size() == size_t(n_correlations_ * n_receivers_)); @@ -288,10 +288,11 @@ void Flagger::computeStats(const std::vector<std::complex<double>>& data) { computeStatisticsWatch_.stop(); } -void Flagger::computeStats(const std::vector<std::complex<float>>& data) { +void Flagger::ComputeStats(const std::vector<std::complex<float>>& data) { computeStatisticsWatch_.start(); - compute_stats_per_antenna(n_receivers_, n_baselines_, n_channels_, - n_correlations_, data, &stats_std_, &stats_sump2_); + detail::ComputeStatsPerAntenna(n_receivers_, n_baselines_, n_channels_, + n_correlations_, data, &stats_std_, + &stats_sump2_); #if defined(DEBUG) assert(stats_std_.size() == size_t(n_correlations_ * n_receivers_)); assert(stats_sump2_.size() == size_t(n_correlations_ * n_receivers_)); @@ -299,26 +300,26 @@ void Flagger::computeStats(const std::vector<std::complex<float>>& data) { computeStatisticsWatch_.stop(); } -void Flagger::assertStatsComputed() const { +void Flagger::AssertStatsComputed() const { if (!(stats_std_.size() == (n_correlations_ * n_receivers_)) || !(stats_sump2_.size() == (n_correlations_ * n_receivers_))) { throw std::runtime_error("Call Flagger::compute_stats() first!"); } } -std::vector<size_t> Flagger::findBadAntennas(float sigma, int maxiters, +std::vector<size_t> Flagger::FindBadAntennas(float sigma, int maxiters, int outlier_threshold) { - assertStatsComputed(); + AssertStatsComputed(); findbadAntennasWatch_.start(); std::vector<size_t> flagged_antennas; - find_bad_antennas_(n_stations_, n_receivers_per_station, n_baselines_, - n_correlations_, sigma, maxiters, stats_sump2_, - &flagged_antennas); - find_bad_antennas_(n_stations_, n_receivers_per_station, n_baselines_, - n_correlations_, sigma, maxiters, stats_std_, - &flagged_antennas); + detail::FindBadAntennas(n_stations_, n_receivers_per_station, n_baselines_, + n_correlations_, sigma, maxiters, stats_sump2_, + &flagged_antennas); + detail::FindBadAntennas(n_stations_, n_receivers_per_station, n_baselines_, + n_correlations_, sigma, maxiters, stats_std_, + &flagged_antennas); // Count the number of times an antenna is flagged std::vector<size_t> flagged_count(n_receivers_); @@ -339,8 +340,8 @@ std::vector<size_t> Flagger::findBadAntennas(float sigma, int maxiters, return flagged_antennas2; } -std::vector<size_t> Flagger::findBadStations(float sigma, int maxiters) { - assertStatsComputed(); +std::vector<size_t> Flagger::FindBadStations(float sigma, int maxiters) { + AssertStatsComputed(); findBadStationsWatch_.start(); @@ -358,8 +359,8 @@ std::vector<size_t> Flagger::findBadStations(float sigma, int maxiters) { } for (size_t i = 0; i < check_sums.size(); ++i) { - std::vector<std::complex<STATS_TYPE>> means = - compute_means(n_stations_, n_receivers_per_station, check_sums[i]); + std::vector<std::complex<STATS_TYPE>> means = detail::ComputeMeans( + n_stations_, n_receivers_per_station, check_sums[i]); for (size_t station = 0; station < n_stations_; ++station) { all_means[0 * check_sums.size() + i][station] = means[station].real(); all_means[1 * check_sums.size() + i][station] = means[station].imag(); diff --git a/antennaflagger/Flagger.h b/antennaflagger/Flagger.h index 1682a0f80..2e2e65452 100644 --- a/antennaflagger/Flagger.h +++ b/antennaflagger/Flagger.h @@ -37,18 +37,18 @@ class Flagger { n_receivers_(n_stations * n_receivers_per_station), n_baselines_(dp3::common::ComputeNBaselines(n_receivers_)) {} - void computeStats(const std::vector<std::complex<double>>& data); - void computeStats(const std::vector<std::complex<float>>& data); + void ComputeStats(const std::vector<std::complex<double>>& data); + void ComputeStats(const std::vector<std::complex<float>>& data); - std::vector<size_t> findBadAntennas(float sigma, int maxiters, + std::vector<size_t> FindBadAntennas(float sigma, int maxiters, int outlier_threshold); - std::vector<size_t> findBadStations(float sigma, int maxiters); + std::vector<size_t> FindBadStations(float sigma, int maxiters); - void reportRuntime() const; + void ReportRuntime() const; private: // Helper function - void assertStatsComputed() const; + void AssertStatsComputed() const; // Constants const size_t n_stations_; diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index c6cee7af6..e43446957 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -202,13 +202,13 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { flaggingTimer_.start(); const std::vector<std::complex<float>> data(buffer_.getData().begin(), buffer_.getData().end()); - flagger_->computeStats(data); + flagger_->ComputeStats(data); // Find outliers both at antenna and station level - const std::vector<size_t> flagged_antennas_1 = flagger_->findBadAntennas( + const std::vector<size_t> flagged_antennas_1 = flagger_->FindBadAntennas( antenna_flagging_sigma_, antenna_flagging_maxiters_, antenna_flagging_outlier_threshold_); - const std::vector<size_t> flagged_antennas_2 = flagger_->findBadStations( + const std::vector<size_t> flagged_antennas_2 = flagger_->FindBadStations( station_flagging_sigma_, station_flagging_maxiters_); // Combine vectors into a single set. -- GitLab From 736f85e4d8dc7ac44ffb39584a15460dcb2472f8 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 10:55:36 +0200 Subject: [PATCH 47/61] Rename DPPP to DP3 --- antennaflagger/Flagger.h | 6 +++--- steps/AntennaFlagger.h | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/antennaflagger/Flagger.h b/antennaflagger/Flagger.h index 2e2e65452..5e24759e2 100644 --- a/antennaflagger/Flagger.h +++ b/antennaflagger/Flagger.h @@ -1,8 +1,8 @@ // Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy) // SPDX-License-Identifier: GPL-3.0-or-later -#ifndef DPPP_ANTENNAFLAGGER_H -#define DPPP_ANTENNAFLAGGER_H +#ifndef DP3_ANTENNAFLAGGER_H +#define DP3_ANTENNAFLAGGER_H #include <vector> #include <complex> @@ -74,4 +74,4 @@ class Flagger { } // namespace antennaflagger } // namespace dp3 -#endif // DPPP_ANTENNAFLAGGER_H +#endif // DP3_ANTENNAFLAGGER_H diff --git a/steps/AntennaFlagger.h b/steps/AntennaFlagger.h index 7bd5048bd..35d1639ae 100644 --- a/steps/AntennaFlagger.h +++ b/steps/AntennaFlagger.h @@ -1,8 +1,8 @@ // Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy) // SPDX-License-Identifier: GPL-3.0-or-later -#ifndef DPPP_ANTENNAFLAGGERSTEP_H -#define DPPP_ANTENNAFLAGGERSTEP_H +#ifndef DP3_ANTENNAFLAGGERSTEP_H +#define DP3_ANTENNAFLAGGERSTEP_H #include <memory> #include <string> @@ -55,4 +55,4 @@ class AntennaFlagger final : public Step { } // namespace steps } // namespace dp3 -#endif // DPPP_ANTENNAFLAGGERSTEP_H +#endif // DP3_ANTENNAFLAGGERSTEP_H -- GitLab From 7df81cab95f85b5f1e2bd36cdffdbb9ca2a78854 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 10:57:36 +0200 Subject: [PATCH 48/61] Remove debug code --- antennaflagger/Flagger.cc | 25 ------------------------- steps/AntennaFlagger.cc | 8 -------- 2 files changed, 33 deletions(-) diff --git a/antennaflagger/Flagger.cc b/antennaflagger/Flagger.cc index e37c39320..9bca21fc5 100644 --- a/antennaflagger/Flagger.cc +++ b/antennaflagger/Flagger.cc @@ -1,9 +1,7 @@ // Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy) // SPDX-License-Identifier: GPL-3.0-or-later -#if defined(DEBUG) #include <cassert> -#endif #include <cstddef> #include <iostream> #include <set> @@ -33,10 +31,8 @@ void ComputeStatsPerBaseline(size_t n_baselines, size_t n_channels, // input and output. This is also different from the Python reference, // and is to avoid strided access later on where the stats are used. -#if defined(DEBUG) size_t n_in = n_baselines * n_channels * n_correlations; assert(data.size() == n_in); -#endif std::vector<OutType> result(n_baselines * n_correlations); size_t n_out = n_correlations * n_baselines; @@ -95,17 +91,13 @@ void ComputeStatsPerAntenna( // - in: data[n_baselines][n_channels][n_correlations] // - out: stats_std[n_correlations][n_baselines] // - out: stats_sump2[n_correlations][n_baselines] -#if defined(DEBUG) assert(data.size() == size_t(n_baselines * n_channels * n_correlations)); -#endif std::vector<std::complex<OutType>> stats_std; std::vector<std::complex<OutType>> stats_sump2; ComputeStatsPerBaseline(n_baselines, n_channels, n_correlations, data, &stats_std, &stats_sump2); -#if defined(DEBUG) assert(stats_std.size() == size_t(n_correlations * n_baselines)); assert(stats_sump2.size() == size_t(n_correlations * n_baselines)); -#endif // Group stats by antenna: // - in: stats_std[n_correlations][n_baselines] @@ -143,10 +135,8 @@ void FindBadAntennas(size_t n_stations, size_t n_receivers_per_station, int n_receivers = n_stations * n_receivers_per_station; // - in: stats[n_correlations][n_receivers] // - out: all_bad_stations[?] -#if defined(DEBUG) assert(n_correlations == static_cast<int>(4)); assert(stats.size() == static_cast<size_t>(n_correlations * n_receivers)); -#endif // Find bad stations by looking at the statistics (per correlation) for all // antennas that belong to one station. The results are stored in a list of @@ -177,22 +167,13 @@ void FindBadAntennas(size_t n_stations, size_t n_receivers_per_station, const size_t n_flags = n_receivers_per_station * 2; std::vector<bool> flags = FindOutliers(sigma, maxiters, check_sum_station); -#if defined(DEBUG) assert(flags.size() == size_t(n_flags)); -#endif // Count number of outliers int n_outliers = 0; for (size_t i = 0; i < n_flags; ++i) { n_outliers += flags[i]; } -#if defined(DEBUG) - if (n_outliers > 0) { - std::cout << "found " << n_outliers - << " outlier antenna(s) for correlation = " << cor - << ", station = " << station << '\n'; - } -#endif // Add the outliers to the list of bad stations const size_t n = bad_stations[station].size(); @@ -281,10 +262,8 @@ void Flagger::ComputeStats(const std::vector<std::complex<double>>& data) { detail::ComputeStatsPerAntenna(n_receivers_, n_baselines_, n_channels_, n_correlations_, data, &stats_std_, &stats_sump2_); -#if defined(DEBUG) assert(stats_std_.size() == size_t(n_correlations_ * n_receivers_)); assert(stats_sump2_.size() == size_t(n_correlations_ * n_receivers_)); -#endif computeStatisticsWatch_.stop(); } @@ -293,10 +272,8 @@ void Flagger::ComputeStats(const std::vector<std::complex<float>>& data) { detail::ComputeStatsPerAntenna(n_receivers_, n_baselines_, n_channels_, n_correlations_, data, &stats_std_, &stats_sump2_); -#if defined(DEBUG) assert(stats_std_.size() == size_t(n_correlations_ * n_receivers_)); assert(stats_sump2_.size() == size_t(n_correlations_ * n_receivers_)); -#endif computeStatisticsWatch_.stop(); } @@ -374,9 +351,7 @@ std::vector<size_t> Flagger::FindBadStations(float sigma, int maxiters) { // Find outliers std::vector<bool> flags = FindOutliers(sigma, maxiters, means); -#if defined(DEBUG) assert(flags.size() == size_t(n_stations_)); -#endif for (size_t station = 0; station < n_stations_; ++station) { station_flags[station] = diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index e43446957..e3547af18 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -217,14 +217,6 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { flagged_antennas.insert(flagged_antennas_1.begin(), flagged_antennas_1.end()); -#if defined(DEBUG) - // Output a list of the flagged antennas - for (int ant : collection) { - std::cout << ant << ","; - } - std::cout << '\n'; -#endif - // Flag all antennas for (int antenna : flagged_antennas) { for (size_t i = 0; i < n_antennas_; i++) { -- GitLab From 278050736248ca6bf50e72622b3416c2a0c37b9f Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 13:55:37 +0200 Subject: [PATCH 49/61] Replace Step::ShPtr by std::shared_ptr --- steps/test/unit/tAntennaFlagger.cc | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/steps/test/unit/tAntennaFlagger.cc b/steps/test/unit/tAntennaFlagger.cc index 15b73c492..7762498fc 100644 --- a/steps/test/unit/tAntennaFlagger.cc +++ b/steps/test/unit/tAntennaFlagger.cc @@ -169,11 +169,14 @@ class TestOutput : public Step { void test1(int ntime, int nbl, int nchan, int ncorr, bool flag) { // Create the steps. TestInput* in = new TestInput(ntime, nbl, nchan, ncorr, flag); - Step::ShPtr step1(in); + std::shared_ptr<TestInput> step1 = + std::make_shared<TestInput>(ntime, nbl, nchan, ncorr, flag); ParameterSet parset; parset.add("selection", "4,10,11"); - Step::ShPtr step2 = std::make_shared<AntennaFlagger>(in, parset, ""); - Step::ShPtr step3 = std::make_shared<TestOutput>(ntime, nbl, nchan, ncorr); + std::shared_ptr<AntennaFlagger> step2 = + std::make_shared<AntennaFlagger>(in, parset, ""); + std::shared_ptr<TestOutput> step3 = + std::make_shared<TestOutput>(ntime, nbl, nchan, ncorr); dp3::steps::test::Execute({step1, step2, step3}); } -- GitLab From 19d6db0cc49fe0903db1fa3b7cc5188866894084 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 14:19:18 +0200 Subject: [PATCH 50/61] Cleanup/fix tests --- common/test/unit/tBaselineUtils.cc | 106 ++++++++++++------------ steps/test/unit/tAntennaFlagger.cc | 127 +++++++++++++++-------------- 2 files changed, 118 insertions(+), 115 deletions(-) diff --git a/common/test/unit/tBaselineUtils.cc b/common/test/unit/tBaselineUtils.cc index 159f183ae..99f902cd8 100644 --- a/common/test/unit/tBaselineUtils.cc +++ b/common/test/unit/tBaselineUtils.cc @@ -5,56 +5,55 @@ BOOST_AUTO_TEST_SUITE(baselineutils) BOOST_AUTO_TEST_CASE(nbaselines) { - for (int nAntennas = 2; nAntennas < 10; nAntennas++) { - int nBaselines = nAntennas * (nAntennas - 1) / 2 + nAntennas; - BOOST_TEST(nBaselines == dp3::common::ComputeNBaselines(nAntennas)); + for (size_t n_antennas = 2; n_antennas < 10; ++n_antennas) { + const size_t n_baselines = n_antennas * (n_antennas - 1) / 2 + n_antennas; + BOOST_TEST(n_baselines == dp3::common::ComputeNBaselines(n_antennas)); } } /* * Helper function to test the row-major baseline utilities */ -void testComputeBaselineRowMajor(int nAntennas) { +void testComputeBaselineRowMajor(size_t n_antennas) { const dp3::common::BaselineOrder order = - dp3::common::BaselineOrder::kRowMajorJOR; - + dp3::common::BaselineOrder::kRowMajor; // Iterate all baselines - int i = 0; - for (int antenna1 = 0; antenna1 < nAntennas; antenna1++) { - for (int antenna2 = antenna1; antenna2 < nAntennas; antenna2++) { + size_t i = 0; + for (size_t antenna1 = 0; antenna1 < n_antennas; ++antenna1) { + for (size_t antenna2 = antenna1; antenna2 < n_antennas; ++antenna2) { // Check reference and computed baseline - std::pair<int, int> baseline_reference = {antenna1, antenna2}; - std::pair<int, int> baseline = - dp3::common::ComputeBaseline(i, nAntennas, order); + const std::pair<size_t, size_t> baseline_reference = {antenna1, antenna2}; + const std::pair<size_t, size_t> baseline = + dp3::common::ComputeBaseline(i, n_antennas, order); BOOST_TEST(baseline_reference.first == baseline.first); BOOST_TEST(baseline_reference.second == baseline.second); // Check baseline indices - int baselineIndex = dp3::common::ComputeBaselineIndex(antenna1, antenna2, - nAntennas, order); - BOOST_TEST(i == baselineIndex); + const size_t baseline_index = dp3::common::ComputeBaselineIndex( + antenna1, antenna2, n_antennas, order); + BOOST_TEST(i == baseline_index); i++; } - // Check baseline indices for antenna1, it should be nAntennas - std::vector<int> baselineIndices = - dp3::common::ComputeBaselineList(antenna1, nAntennas, order); - BOOST_TEST(int(baselineIndices.size()) == nAntennas); + // Check baseline indices for antenna1, it should be n_antennas + const std::vector<size_t> baselineIndices = + dp3::common::ComputeBaselineList(antenna1, n_antennas, order); + BOOST_TEST(baselineIndices.size() == n_antennas); // Use the baseline indices to lookup the baseline and check whether - // antenna1 occurs exactly nAntennas + 1 times (correlations + 1 // + // antenna1 occurs exactly n_antennas + 1 times (correlations + 1 // autocorrelation). - int count = 0; - for (int antenna2 = 0; antenna2 < nAntennas; antenna2++) { - const int baselineIndex = dp3::common::ComputeBaselineIndex( - antenna1, antenna2, nAntennas, order); - BOOST_TEST(baselineIndices[antenna2] == baselineIndex); - std::pair<int, int> baseline = - dp3::common::ComputeBaseline(baselineIndex, nAntennas, order); + size_t count = 0; + for (size_t antenna2 = 0; antenna2 < n_antennas; ++antenna2) { + const size_t baseline_index = dp3::common::ComputeBaselineIndex( + antenna1, antenna2, n_antennas, order); + BOOST_TEST(baselineIndices[antenna2] == baseline_index); + const std::pair<size_t, size_t> baseline = + dp3::common::ComputeBaseline(baseline_index, n_antennas, order); count += baseline.first == antenna1; count += baseline.second == antenna1; } - BOOST_TEST(count == nAntennas + 1); + BOOST_TEST(count == n_antennas + 1); } } @@ -66,47 +65,46 @@ BOOST_AUTO_TEST_CASE(rowmajor) { /* * Helper function to test the column-major baseline utilities */ -void testComputeBaselineColumnMajor(int nAntennas) { +void testComputeBaselineColumnMajor(size_t n_antennas) { const dp3::common::BaselineOrder order = - dp3::common::BaselineOrder::COLUMN_MAJOR; - + dp3::common::BaselineOrder::kColumnMajor; // Iterate all baselines - int i = 0; - for (int antenna2 = 0; antenna2 < nAntennas; antenna2++) { - for (int antenna1 = 0; antenna1 <= antenna2; antenna1++) { + size_t i = 0; + for (size_t antenna2 = 0; antenna2 < n_antennas; ++antenna2) { + for (size_t antenna1 = 0; antenna1 <= antenna2; ++antenna1) { // Check reference and computed baseline - std::pair<int, int> baseline_reference = {antenna1, antenna2}; - std::pair<int, int> baseline = - dp3::common::ComputeBaseline(i, nAntennas, order); + std::pair<size_t, size_t> baseline_reference = {antenna1, antenna2}; + std::pair<size_t, size_t> baseline = + dp3::common::ComputeBaseline(i, n_antennas, order); BOOST_TEST(baseline_reference.first == baseline.first); BOOST_TEST(baseline_reference.second == baseline.second); // Check baseline indices - int baselineIndex = dp3::common::ComputeBaselineIndex(antenna1, antenna2, - nAntennas, order); - BOOST_TEST(i == baselineIndex); - i++; + size_t baseline_index = dp3::common::ComputeBaselineIndex( + antenna1, antenna2, n_antennas, order); + BOOST_TEST(i == baseline_index); + ++i; } - // Check baseline indices for antenna2, it should be nAntennas - std::vector<int> baselineIndices = - dp3::common::ComputeBaselineList(antenna2, nAntennas, order); - BOOST_TEST(int(baselineIndices.size()) == nAntennas); + // Check baseline indices for antenna2, it should be n_antennas + const std::vector<size_t> baseline_indices = + dp3::common::ComputeBaselineList(antenna2, n_antennas, order); + BOOST_TEST(baseline_indices.size() == n_antennas); // Use the baseline indices to lookup the baseline and check whether - // antenna1 occurs exactly nAntennas + 1 times (correlations + 1 // + // antenna1 occurs exactly n_antennas + 1 times (correlations + 1 // autocorrelation). - int count = 0; - for (int antenna1 = 0; antenna1 < nAntennas; antenna1++) { - const int baselineIndex = dp3::common::ComputeBaselineIndex( - antenna1, antenna2, nAntennas, order); - BOOST_TEST(baselineIndices[antenna1] == baselineIndex); - std::pair<int, int> baseline = - dp3::common::ComputeBaseline(baselineIndex, nAntennas, order); + size_t count = 0; + for (size_t antenna1 = 0; antenna1 < n_antennas; ++antenna1) { + const size_t baseline_index = dp3::common::ComputeBaselineIndex( + antenna1, antenna2, n_antennas, order); + BOOST_TEST(baseline_indices[antenna1] == baseline_index); + const std::pair<size_t, size_t> baseline = + dp3::common::ComputeBaseline(baseline_index, n_antennas, order); count += baseline.first == antenna2; count += baseline.second == antenna2; } - BOOST_TEST(count == nAntennas + 1); + BOOST_TEST(count == n_antennas + 1); } } diff --git a/steps/test/unit/tAntennaFlagger.cc b/steps/test/unit/tAntennaFlagger.cc index 7762498fc..bfb04abac 100644 --- a/steps/test/unit/tAntennaFlagger.cc +++ b/steps/test/unit/tAntennaFlagger.cc @@ -32,28 +32,30 @@ BOOST_AUTO_TEST_SUITE(antennaflagger) // It can be used with different nr of times, channels, etc. class TestInput final : public InputStep { public: - TestInput(int ntime, int nbl, int nchan, int ncorr, bool flag) - : count_(0), - n_times_(ntime), - n_bl_(nbl), - n_chan_(nchan), - n_corr__(ncorr), + TestInput(int n_time, int n_baselines, int n_channels, int n_correlations, + bool flag) + : n_time_processed_(0), + n_time_(n_time), + n_baselines_(n_baselines), + n_channels_(n_channels), + n_correlations_(n_correlations), flag_(flag) { // Define start time 0.5 (= 3 - 0.5*5) and time interval 5. - info().init(ncorr, 0, nchan, ntime, 0.5, 5., string(), string()); + info().init(n_correlations, 0, n_channels, n_time, 0.5, 5., string(), + string()); // Fill the baselines; use 12 stations with 48 antennas each. - constexpr int nreceivers = 12 * 48; - std::vector<int> ant1(nbl); - std::vector<int> ant2(nbl); + constexpr int n_receivers = 12 * 48; + std::vector<int> ant1(n_baselines); + std::vector<int> ant2(n_baselines); int st1 = 0; int st2 = 0; - for (int i = 0; i < nbl; ++i) { + for (int i = 0; i < n_baselines; ++i) { ant1[i] = st1; ant2[i] = st2; - if (++st2 == nreceivers) { + if (++st2 == n_receivers) { st2 = 0; - if (++st1 == nreceivers) { + if (++st1 == n_receivers) { st1 = 0; } } @@ -61,45 +63,48 @@ class TestInput final : public InputStep { // Atenna names and positions are not used by the AntennaFlagger, // but the arrays need to have the correct dimensions. - std::vector<std::string> antNames(nreceivers); - std::vector<casacore::MPosition> antPos(nreceivers); - std::vector<double> antDiam(nreceivers, 70.); + std::vector<std::string> antNames(n_receivers); + std::vector<casacore::MPosition> antPos(n_receivers); + std::vector<double> antDiam(n_receivers, 70.); info().set(antNames, antDiam, antPos, ant1, ant2); } private: bool process(const DPBuffer&) override { // Stop when all times are done. - if (count_ == n_times_) { + if (n_time_processed_ == n_time_) { return false; } - casacore::Cube<casacore::Complex> data(n_corr__, n_chan_, n_bl_); + casacore::Cube<casacore::Complex> data(n_correlations_, n_channels_, + n_baselines_); for (size_t i = 0; i < data.size(); ++i) { - data.data()[i] = casacore::Complex(i + count_ * 10, i - 10 + count_ * 6); + data.data()[i] = casacore::Complex(i + n_time_processed_ * 10, + i - 10 + n_time_processed_ * 6); } - casacore::Matrix<double> uvw(3, n_bl_); - for (int i = 0; i < n_bl_; ++i) { - uvw(0, i) = 1 + count_ + i; - uvw(1, i) = 2 + count_ + i; - uvw(2, i) = 3 + count_ + i; + casacore::Matrix<double> uvw(3, n_baselines_); + for (size_t i = 0; i < n_baselines_; ++i) { + uvw(0, i) = 1 + n_time_processed_ + i; + uvw(1, i) = 2 + n_time_processed_ + i; + uvw(2, i) = 3 + n_time_processed_ + i; } - DPBuffer buf; - buf.setTime(count_ * 5 + 3); // same interval as in updateAveragInfo - buf.setData(data); - buf.setUVW(uvw); + DPBuffer buffer; + buffer.setTime(n_time_processed_ * 5 + + 3); // same interval as in updateAveragInfo + buffer.setData(data); + buffer.setUVW(uvw); casacore::Cube<float> weights(data.shape()); weights = 1.0; - buf.setWeights(weights); + buffer.setWeights(weights); casacore::Cube<bool> flags(data.shape()); flags = flag_; - buf.setFlags(flags); + buffer.setFlags(flags); // The fullRes flags are a copy of the XX flags, but differently shaped. // They are not averaged, thus only 1 time per row. - casacore::Cube<bool> fullResFlags(n_chan_, 1, n_bl_); + casacore::Cube<bool> fullResFlags(n_channels_, 1, n_baselines_); fullResFlags = flag_; - buf.setFullResFlags(fullResFlags); - getNextStep()->process(buf); - ++count_; + buffer.setFullResFlags(fullResFlags); + getNextStep()->process(buffer); + ++n_time_processed_; return true; } @@ -107,40 +112,40 @@ class TestInput final : public InputStep { void show(std::ostream&) const override{}; void updateInfo(const DPInfo&) override {} - int count_; - int n_times_; - int n_bl_; - int n_chan_; - int n_corr__; + size_t n_time_processed_; + size_t n_time_; + size_t n_baselines_; + size_t n_channels_; + size_t n_correlations_; bool flag_; }; // Class to check result of flagged, unaveraged TestInput run by test1. class TestOutput : public Step { public: - TestOutput(int ntime, int nbl, int nchan, int ncorr) - : count_(0), - n_times_(ntime), - n_bl_(nbl), - n_chan_(nchan), - n_corr_(ncorr) {} + TestOutput(int n_time, int n_baselines, int n_channels, int n_correlations) + : n_time_processed_(0), + n_time_(n_time), + n_baselines_(n_baselines), + n_channels_(n_channels), + n_correlations_(n_correlations) {} private: - virtual bool process(const DPBuffer& buf) { - casacore::Cube<bool> result(n_corr_, n_chan_, n_bl_); - for (int i = 0; i < n_bl_; ++i) { - casacore::Cube<bool> test = buf.getFlags(); + virtual bool process(const DPBuffer& buffer) { + casacore::Cube<bool> result(n_correlations_, n_channels_, n_baselines_); + for (int i = 0; i < n_baselines_; ++i) { + casacore::Cube<bool> test = buffer.getFlags(); if (i == 4 || i == 10 || i == 11) { - for (int j = 0; j < n_chan_; ++j) { - for (int k = 0; k < n_corr_; ++k) { + for (int j = 0; j < n_channels_; ++j) { + for (int k = 0; k < n_correlations_; ++k) { result(k, j, i) = true; } } } } - BOOST_CHECK(allEQ(buf.getFlags(), result)); - count_++; + BOOST_CHECK(allEQ(buffer.getFlags(), result)); + n_time_processed_++; return true; } @@ -148,19 +153,19 @@ class TestOutput : public Step { virtual void show(std::ostream&) const {} virtual void updateInfo(const DPInfo& infoIn) { info() = infoIn; - BOOST_CHECK_EQUAL(static_cast<int>(infoIn.origNChan()), n_chan_); - BOOST_CHECK_EQUAL(static_cast<int>(infoIn.nchan()), n_chan_); - BOOST_CHECK_EQUAL(static_cast<int>(infoIn.ntime()), n_times_); + BOOST_CHECK_EQUAL(static_cast<int>(infoIn.origNChan()), n_channels_); + BOOST_CHECK_EQUAL(static_cast<int>(infoIn.nchan()), n_channels_); + BOOST_CHECK_EQUAL(static_cast<int>(infoIn.ntime()), n_time_); BOOST_CHECK_EQUAL(infoIn.timeInterval(), 5); BOOST_CHECK_EQUAL(static_cast<int>(infoIn.nchanAvg()), 1); BOOST_CHECK_EQUAL(static_cast<int>(infoIn.ntimeAvg()), 1); } - int count_; - int n_times_; - int n_bl_; - int n_chan_; - int n_corr_; + int n_time_processed_; + int n_time_; + int n_baselines_; + int n_channels_; + int n_correlations_; bool flag_; bool clear_; }; -- GitLab From d3d20ee14b6b04f39fd574f413e7b94785c3ee82 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 14:21:28 +0200 Subject: [PATCH 51/61] Update tBaselineUtils --- common/test/unit/tBaselineUtils.cc | 32 +++++++++++++++--------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/common/test/unit/tBaselineUtils.cc b/common/test/unit/tBaselineUtils.cc index 99f902cd8..b10f73aa3 100644 --- a/common/test/unit/tBaselineUtils.cc +++ b/common/test/unit/tBaselineUtils.cc @@ -18,21 +18,21 @@ void testComputeBaselineRowMajor(size_t n_antennas) { const dp3::common::BaselineOrder order = dp3::common::BaselineOrder::kRowMajor; // Iterate all baselines - size_t i = 0; + size_t baseline_count = 0; for (size_t antenna1 = 0; antenna1 < n_antennas; ++antenna1) { for (size_t antenna2 = antenna1; antenna2 < n_antennas; ++antenna2) { // Check reference and computed baseline const std::pair<size_t, size_t> baseline_reference = {antenna1, antenna2}; const std::pair<size_t, size_t> baseline = - dp3::common::ComputeBaseline(i, n_antennas, order); + dp3::common::ComputeBaseline(baseline_count, n_antennas, order); BOOST_TEST(baseline_reference.first == baseline.first); BOOST_TEST(baseline_reference.second == baseline.second); // Check baseline indices const size_t baseline_index = dp3::common::ComputeBaselineIndex( antenna1, antenna2, n_antennas, order); - BOOST_TEST(i == baseline_index); - i++; + BOOST_TEST(baseline_count == baseline_index); + ++baseline_count; } // Check baseline indices for antenna1, it should be n_antennas @@ -43,17 +43,17 @@ void testComputeBaselineRowMajor(size_t n_antennas) { // Use the baseline indices to lookup the baseline and check whether // antenna1 occurs exactly n_antennas + 1 times (correlations + 1 // autocorrelation). - size_t count = 0; + size_t antenna_count = 0; for (size_t antenna2 = 0; antenna2 < n_antennas; ++antenna2) { const size_t baseline_index = dp3::common::ComputeBaselineIndex( antenna1, antenna2, n_antennas, order); BOOST_TEST(baselineIndices[antenna2] == baseline_index); const std::pair<size_t, size_t> baseline = dp3::common::ComputeBaseline(baseline_index, n_antennas, order); - count += baseline.first == antenna1; - count += baseline.second == antenna1; + antenna_count += baseline.first == antenna1; + antenna_count += baseline.second == antenna1; } - BOOST_TEST(count == n_antennas + 1); + BOOST_TEST(antenna_count == n_antennas + 1); } } @@ -69,21 +69,21 @@ void testComputeBaselineColumnMajor(size_t n_antennas) { const dp3::common::BaselineOrder order = dp3::common::BaselineOrder::kColumnMajor; // Iterate all baselines - size_t i = 0; + size_t baseline_count = 0; for (size_t antenna2 = 0; antenna2 < n_antennas; ++antenna2) { for (size_t antenna1 = 0; antenna1 <= antenna2; ++antenna1) { // Check reference and computed baseline std::pair<size_t, size_t> baseline_reference = {antenna1, antenna2}; std::pair<size_t, size_t> baseline = - dp3::common::ComputeBaseline(i, n_antennas, order); + dp3::common::ComputeBaseline(baseline_count, n_antennas, order); BOOST_TEST(baseline_reference.first == baseline.first); BOOST_TEST(baseline_reference.second == baseline.second); // Check baseline indices size_t baseline_index = dp3::common::ComputeBaselineIndex( antenna1, antenna2, n_antennas, order); - BOOST_TEST(i == baseline_index); - ++i; + BOOST_TEST(baseline_count == baseline_index); + ++baseline_count; } // Check baseline indices for antenna2, it should be n_antennas @@ -94,17 +94,17 @@ void testComputeBaselineColumnMajor(size_t n_antennas) { // Use the baseline indices to lookup the baseline and check whether // antenna1 occurs exactly n_antennas + 1 times (correlations + 1 // autocorrelation). - size_t count = 0; + size_t antenna_count = 0; for (size_t antenna1 = 0; antenna1 < n_antennas; ++antenna1) { const size_t baseline_index = dp3::common::ComputeBaselineIndex( antenna1, antenna2, n_antennas, order); BOOST_TEST(baseline_indices[antenna1] == baseline_index); const std::pair<size_t, size_t> baseline = dp3::common::ComputeBaseline(baseline_index, n_antennas, order); - count += baseline.first == antenna2; - count += baseline.second == antenna2; + antenna_count += baseline.first == antenna2; + antenna_count += baseline.second == antenna2; } - BOOST_TEST(count == n_antennas + 1); + BOOST_TEST(antenna_count == n_antennas + 1); } } -- GitLab From 0f9dbe4eb1ce87856c873c153eb357d419e5d164 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 14:24:01 +0200 Subject: [PATCH 52/61] Remove unneeded inline --- common/baseline_indices.h | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/common/baseline_indices.h b/common/baseline_indices.h index 763d6e232..132cf0158 100644 --- a/common/baseline_indices.h +++ b/common/baseline_indices.h @@ -91,9 +91,9 @@ inline size_t ComputeBaselineIndex(size_t antenna_a, size_t antenna_b, * Computes the antenna pair (baseline) given the baseline index. It assumes * that auto-correlations of antennas are included. */ -inline std::pair<size_t, size_t> ComputeBaseline(size_t baseline_index, - size_t n_antennas, - BaselineOrder order) { +std::pair<size_t, size_t> ComputeBaseline(size_t baseline_index, + size_t n_antennas, + BaselineOrder order) { if (order == BaselineOrder::kRowMajor) { const size_t n_baselines = ComputeNBaselines(n_antennas); const size_t antenna = @@ -116,9 +116,8 @@ inline std::pair<size_t, size_t> ComputeBaseline(size_t baseline_index, * Computes the list of baselines for a given antenna. It assumes that * auto-correlations are included. */ -inline std::vector<size_t> ComputeBaselineList(size_t antenna, - size_t n_antennas, - BaselineOrder order) { +std::vector<size_t> ComputeBaselineList(size_t antenna, size_t n_antennas, + BaselineOrder order) { std::vector<size_t> indices; indices.reserve(n_antennas); -- GitLab From 31195b42e4cb45e948cf4a85b991f7a698e1cdec Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 16:33:32 +0200 Subject: [PATCH 53/61] Revert inline to workaround compilation errors --- common/baseline_indices.h | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/common/baseline_indices.h b/common/baseline_indices.h index 132cf0158..763d6e232 100644 --- a/common/baseline_indices.h +++ b/common/baseline_indices.h @@ -91,9 +91,9 @@ inline size_t ComputeBaselineIndex(size_t antenna_a, size_t antenna_b, * Computes the antenna pair (baseline) given the baseline index. It assumes * that auto-correlations of antennas are included. */ -std::pair<size_t, size_t> ComputeBaseline(size_t baseline_index, - size_t n_antennas, - BaselineOrder order) { +inline std::pair<size_t, size_t> ComputeBaseline(size_t baseline_index, + size_t n_antennas, + BaselineOrder order) { if (order == BaselineOrder::kRowMajor) { const size_t n_baselines = ComputeNBaselines(n_antennas); const size_t antenna = @@ -116,8 +116,9 @@ std::pair<size_t, size_t> ComputeBaseline(size_t baseline_index, * Computes the list of baselines for a given antenna. It assumes that * auto-correlations are included. */ -std::vector<size_t> ComputeBaselineList(size_t antenna, size_t n_antennas, - BaselineOrder order) { +inline std::vector<size_t> ComputeBaselineList(size_t antenna, + size_t n_antennas, + BaselineOrder order) { std::vector<size_t> indices; indices.reserve(n_antennas); -- GitLab From 17b50b7e547c8fa35f276d4eca9905b5eb6e8f6b Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 16:56:07 +0200 Subject: [PATCH 54/61] Replace STATS_TYPE definition by StatsType alias --- antennaflagger/Flagger.cc | 14 +++++++------- antennaflagger/Flagger.h | 6 +++--- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/antennaflagger/Flagger.cc b/antennaflagger/Flagger.cc index 9bca21fc5..9bf8ca261 100644 --- a/antennaflagger/Flagger.cc +++ b/antennaflagger/Flagger.cc @@ -324,19 +324,19 @@ std::vector<size_t> Flagger::FindBadStations(float sigma, int maxiters) { // Compute means for all the check sums seperately for // the real and imaginary part of the check sum. - // - in: STATS_TYPE2 check_sums[n_check_sums][n_receivers] - // - out: STATS_TYPE all_means[check_sums.size() * 2][n_stations] - std::array<const std::complex<STATS_TYPE>*, 4> check_sums{ + // - in: StatsType2 check_sums[n_check_sums][n_receivers] + // - out: StatsType all_means[check_sums.size() * 2][n_stations] + std::array<const std::complex<StatsType>*, 4> check_sums{ &stats_std_[0 * n_receivers_], &stats_std_[3 * n_receivers_], &stats_sump2_[0 * n_receivers_], &stats_sump2_[3 * n_receivers_]}; - std::vector<std::vector<STATS_TYPE>> all_means(check_sums.size() * 2); + std::vector<std::vector<StatsType>> all_means(check_sums.size() * 2); - for (std::vector<STATS_TYPE>& means : all_means) { + for (std::vector<StatsType>& means : all_means) { means.resize(n_stations_); } for (size_t i = 0; i < check_sums.size(); ++i) { - std::vector<std::complex<STATS_TYPE>> means = detail::ComputeMeans( + std::vector<std::complex<StatsType>> means = detail::ComputeMeans( n_stations_, n_receivers_per_station, check_sums[i]); for (size_t station = 0; station < n_stations_; ++station) { all_means[0 * check_sums.size() + i][station] = means[station].real(); @@ -347,7 +347,7 @@ std::vector<size_t> Flagger::FindBadStations(float sigma, int maxiters) { // Flag all stations for which all of the means are an outlier std::vector<bool> station_flags(n_stations_, {false}); for (size_t i = 0; i < all_means.size(); ++i) { - const std::vector<STATS_TYPE>& means = all_means[i]; + const std::vector<StatsType>& means = all_means[i]; // Find outliers std::vector<bool> flags = FindOutliers(sigma, maxiters, means); diff --git a/antennaflagger/Flagger.h b/antennaflagger/Flagger.h index 5e24759e2..7360c827b 100644 --- a/antennaflagger/Flagger.h +++ b/antennaflagger/Flagger.h @@ -66,9 +66,9 @@ class Flagger { common::NSTimer findBadStationsWatch_; // Data -#define STATS_TYPE float - std::vector<std::complex<STATS_TYPE>> stats_std_; - std::vector<std::complex<STATS_TYPE>> stats_sump2_; + using StatsType = float; + std::vector<std::complex<StatsType>> stats_std_; + std::vector<std::complex<StatsType>> stats_sump2_; }; } // namespace antennaflagger -- GitLab From 89ceaceeabc23a22c2b4601d3ee63c6b5a54f9c5 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Fri, 19 Aug 2022 17:14:38 +0200 Subject: [PATCH 55/61] Update Flagger and its documentation --- antennaflagger/Flagger.cc | 8 ++++---- antennaflagger/Flagger.h | 3 +-- docs/schemas/AntennaFlagger.yml | 30 ++++++++++++++++-------------- steps/AntennaFlagger.cc | 5 +---- 4 files changed, 22 insertions(+), 24 deletions(-) diff --git a/antennaflagger/Flagger.cc b/antennaflagger/Flagger.cc index 9bf8ca261..63ee20ccd 100644 --- a/antennaflagger/Flagger.cc +++ b/antennaflagger/Flagger.cc @@ -284,8 +284,7 @@ void Flagger::AssertStatsComputed() const { } } -std::vector<size_t> Flagger::FindBadAntennas(float sigma, int maxiters, - int outlier_threshold) { +std::vector<size_t> Flagger::FindBadAntennas(float sigma, int maxiters) { AssertStatsComputed(); findbadAntennasWatch_.start(); @@ -304,10 +303,11 @@ std::vector<size_t> Flagger::FindBadAntennas(float sigma, int maxiters, flagged_count[antenna]++; } - // Make a list of antennas that are flagged multiple times + // Make a list of antennas that are flagged using both of the statistics that + // are used above. std::vector<size_t> flagged_antennas2; for (size_t antenna = 0; antenna < n_receivers_; ++antenna) { - if (flagged_count[antenna] > size_t(outlier_threshold)) { + if (flagged_count[antenna] > size_t(2)) { flagged_antennas2.push_back(antenna); } } diff --git a/antennaflagger/Flagger.h b/antennaflagger/Flagger.h index 7360c827b..fdc946d52 100644 --- a/antennaflagger/Flagger.h +++ b/antennaflagger/Flagger.h @@ -40,8 +40,7 @@ class Flagger { void ComputeStats(const std::vector<std::complex<double>>& data); void ComputeStats(const std::vector<std::complex<float>>& data); - std::vector<size_t> FindBadAntennas(float sigma, int maxiters, - int outlier_threshold); + std::vector<size_t> FindBadAntennas(float sigma, int maxiters); std::vector<size_t> FindBadStations(float sigma, int maxiters); void ReportRuntime() const; diff --git a/docs/schemas/AntennaFlagger.yml b/docs/schemas/AntennaFlagger.yml index 27f82736d..660a25eec 100644 --- a/docs/schemas/AntennaFlagger.yml +++ b/docs/schemas/AntennaFlagger.yml @@ -1,8 +1,9 @@ description: >- - Flag all baselines containing bad antennas. The antennas can be specified by - their number (see `selection`), and/or can be determined automatically based - on (per time slot) statistics when `detect_outliers=True`. This step is mainly - intended for processing of AARTFAAC data. + Flag all baselines containing containing manually selected or automatically + determined antennas. The antennas can be specified by their number (see + `selection`), and/or can be determined automatically based on (per time slot) + statistics when `detect_outliers=True`. This step is mainly intended for + processing of AARTFAAC data. inputs: step_name: type: string @@ -19,32 +20,33 @@ inputs: The antenna numbers (not names) to flag, currently the following selection are implemented: Single antennas (e.g. 0,3,15,4), antenna ranges (e.g. 1~15), autocorrelations for a selection (e.g. 0~3&&& flags all autocorrelations of antennas 0, 1, 2, 3). - Multiple selections can be separated with a semicolon (;) `.` + Multiple selections can be separated with a semicolon (;). Using this + selection is only relevant for AARTFAAC data, where antennas in a station are + processed as seperate receivers `.` detect_outliers: + default: false type: boolean doc: >- If true, outlier detection on antenna and station level is enabled and affected baselines are flagged. Flagging is applied to all antennas, regardless of whether they are in `selection` or not `.` antenna_flagging_sigma: + default: 3 type: float doc: >- - Flag an atenna when its data exceeds this threshold `.` + Flag an antenna when the statistical properties of its data exceeds this threshold. antenna_flagging_maxiters: + default: 5 type: int doc: >- Maximum number of iterations to find outliers in the data `.` - antenna_flagging_outlier_threshold: - type: int - doc: >- - An antenna is flagged when the sum of squares and/or standard deviation of - its data exceeds a threshold. This parameter controls whether either one of - these thresholds or both must exceeded before the antenna is considered to - be an outlier `.` station_flagging_sigma: + default: 2.5 type: float doc: >- - Flag all antennas for this station when its data exceeds this threshold `.` + Flag all antennas for this station when the statistical properties of its + data exceeds this threshold `.` station_flagging_maxiters: + default: 5 type: int doc: >- Maximum number of iterations to find outliers in the data `.` diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index e3547af18..7749104e9 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -147,8 +147,6 @@ AntennaFlagger::AntennaFlagger(InputStep* input, parset.getFloat(prefix + "antenna_flagging_sigma", 3); antenna_flagging_maxiters_ = parset.getInt(prefix + "antenna_flagging_maxiters", 5); - antenna_flagging_outlier_threshold_ = - parset.getInt(prefix + "antenna_flagging_outlier_threshold", 2); station_flagging_sigma_ = parset.getFloat(prefix + "station_flagging_sigma", 2.5); station_flagging_maxiters_ = @@ -206,8 +204,7 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { // Find outliers both at antenna and station level const std::vector<size_t> flagged_antennas_1 = flagger_->FindBadAntennas( - antenna_flagging_sigma_, antenna_flagging_maxiters_, - antenna_flagging_outlier_threshold_); + antenna_flagging_sigma_, antenna_flagging_maxiters_); const std::vector<size_t> flagged_antennas_2 = flagger_->FindBadStations( station_flagging_sigma_, station_flagging_maxiters_); -- GitLab From 997597041923917079a203389a97ac7960e452fd Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Mon, 22 Aug 2022 08:47:20 +0200 Subject: [PATCH 56/61] Replace IsValid by std::isfinite --- antennaflagger/statistics.h | 14 +------------- 1 file changed, 1 insertion(+), 13 deletions(-) diff --git a/antennaflagger/statistics.h b/antennaflagger/statistics.h index f95376208..3e8b50794 100644 --- a/antennaflagger/statistics.h +++ b/antennaflagger/statistics.h @@ -5,24 +5,12 @@ namespace detail { -template <typename T> -inline bool IsValid(T x) { - bool is_nan = x != x; - bool is_finite = x != 2 * x; - return !is_nan && is_finite; -} - -template <typename T> -inline bool IsValid(const std::complex<T>& x) { - return IsValid(x.real()) && IsValid(x.imag()); -} - template <typename T> void FlagInvalid(const std::vector<T>& data, std::vector<bool>& flags) { assert(data.size() == flags.size()); for (size_t i = 0; i < data.size(); ++i) { - if (!IsValid(data[i])) { + if (!std::isfinite(data[i])) { flags[i] = true; } } -- GitLab From b3123220d19c548777b60e5bb94b4bf13e0b5f47 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Tue, 23 Aug 2022 15:45:51 +0200 Subject: [PATCH 57/61] Various style-related changes --- steps/AntennaFlagger.cc | 57 +++++++++++++++--------------- steps/AntennaFlagger.h | 10 +++--- steps/test/unit/tAntennaFlagger.cc | 16 ++++----- 3 files changed, 42 insertions(+), 41 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index 7749104e9..5323c9cbf 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -108,13 +108,13 @@ std::vector<std::pair<int, int>> ConvertSelection( return selection; } -xt::xtensor<bool, 2> Convert(const std::string& selectionString, - size_t nAntennas) { - assert(nAntennas != 0); - xt::xtensor<bool, 2> bl({nAntennas, nAntennas}, false); +xt::xtensor<bool, 2> Convert(const std::string& selection_string, + size_t n_antennas) { + assert(n_antennas != 0); + xt::xtensor<bool, 2> bl({n_antennas, n_antennas}, false); std::vector<std::pair<int, int>> selection = - ConvertSelection(selectionString, static_cast<int>(nAntennas)); + ConvertSelection(selection_string, static_cast<int>(n_antennas)); for (std::pair<int, int> p : selection) { bl(p.first, p.second) = true; @@ -132,17 +132,18 @@ void AntennaFlagger::finish() { getNextStep()->finish(); } void AntennaFlagger::show(std::ostream& ostream) const { ostream << "AntennaFlagger " << name_; - ostream << "\n selection: " << selectionString_; + ostream << "\n selection: " << selection_string_; } AntennaFlagger::AntennaFlagger(InputStep* input, const common::ParameterSet& parset, const string& prefix) { - initializationTimer_.start(); + initialization_timer_.start(); + // Parse arguments - std::string selectionString = + const std::string selection_string = parset.getString(prefix + "selection", std::string()); - doDetectOutliers_ = parset.getBool(prefix + "detect_outliers", false); + do_detect_outliers_ = parset.getBool(prefix + "detect_outliers", false); antenna_flagging_sigma_ = parset.getFloat(prefix + "antenna_flagging_sigma", 3); antenna_flagging_maxiters_ = @@ -156,11 +157,11 @@ AntennaFlagger::AntennaFlagger(InputStep* input, n_antennas_ = input->getInfo().nantenna(); n_correlations_ = input->getInfo().ncorr(); n_channels_ = input->getInfo().nchan(); - selectionString_ = selectionString; - selection_ = Convert(selectionString, n_antennas_); + selection_string_ = selection_string; + selection_ = Convert(selection_string, n_antennas_); name_ = prefix; - if (doDetectOutliers_) { + if (do_detect_outliers_) { std::string antennaSet(input->getInfo().antennaSet()); std::vector<std::string> antennaNames(input->getInfo().antennaNames()); size_t n_receivers_per_station = 0; @@ -181,11 +182,11 @@ AntennaFlagger::AntennaFlagger(InputStep* input, n_stations, n_receivers_per_station, n_channels_, n_correlations_); } - initializationTimer_.stop(); + initialization_timer_.stop(); } bool AntennaFlagger::process(const base::DPBuffer& buffer) { - computationTimer_.start(); + computation_timer_.start(); buffer_.copy(buffer); @@ -195,9 +196,9 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { /* * Do additional flagging on the input data based on outlier statistics. */ - if (doDetectOutliers_) { - computationTimer_.stop(); - flaggingTimer_.start(); + if (do_detect_outliers_) { + computation_timer_.stop(); + flagging_timer_.start(); const std::vector<std::complex<float>> data(buffer_.getData().begin(), buffer_.getData().end()); flagger_->ComputeStats(data); @@ -222,8 +223,8 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { } } - flaggingTimer_.stop(); - computationTimer_.start(); + flagging_timer_.stop(); + computation_timer_.start(); } for (size_t correlation = 0; correlation < flags.nrow(); correlation++) { @@ -239,26 +240,26 @@ bool AntennaFlagger::process(const base::DPBuffer& buffer) { buffer_.setFlags(flags); - computationTimer_.stop(); + computation_timer_.stop(); getNextStep()->process(buffer_); return true; } void AntennaFlagger::showTimings(std::ostream& ostream, double duration) const { - double initTime = initializationTimer_.getElapsed(); - double flagTime = flaggingTimer_.getElapsed(); - double computeTime = computationTimer_.getElapsed(); - double totTime = initTime + flagTime + computeTime; + double initializion_time = initialization_timer_.getElapsed(); + double flagging_time = flagging_timer_.getElapsed(); + double computation_time = computation_timer_.getElapsed(); + double total_time = initializion_time + flagging_time + computation_time; ostream << " "; - base::FlagCounter::showPerc1(ostream, totTime, duration); + base::FlagCounter::showPerc1(ostream, total_time, duration); ostream << " AntennaFlagger " << name_ << "\n "; - base::FlagCounter::showPerc1(ostream, initTime, totTime); + base::FlagCounter::showPerc1(ostream, initializion_time, total_time); ostream << " of it spent in initialization.\n "; - base::FlagCounter::showPerc1(ostream, flagTime, totTime); + base::FlagCounter::showPerc1(ostream, flagging_time, total_time); ostream << " of it spent in flagging.\n "; - base::FlagCounter::showPerc1(ostream, computeTime, totTime); + base::FlagCounter::showPerc1(ostream, computation_time, total_time); ostream << " of it spent in setting flags.\n"; } diff --git a/steps/AntennaFlagger.h b/steps/AntennaFlagger.h index 35d1639ae..a072c97d1 100644 --- a/steps/AntennaFlagger.h +++ b/steps/AntennaFlagger.h @@ -35,8 +35,8 @@ class AntennaFlagger final : public Step { xt::xtensor<bool, 2> selection_; std::string name_; - std::string selectionString_; - bool doDetectOutliers_; + std::string selection_string_; + bool do_detect_outliers_; std::unique_ptr<dp3::antennaflagger::Flagger> flagger_; // AartfaacFlagger parameters @@ -47,9 +47,9 @@ class AntennaFlagger final : public Step { size_t station_flagging_maxiters_; // TImers - common::NSTimer initializationTimer_; - common::NSTimer computationTimer_; - common::NSTimer flaggingTimer_; + common::NSTimer initialization_timer_; + common::NSTimer computation_timer_; + common::NSTimer flagging_timer_; }; } // namespace steps diff --git a/steps/test/unit/tAntennaFlagger.cc b/steps/test/unit/tAntennaFlagger.cc index bfb04abac..2fa5ca973 100644 --- a/steps/test/unit/tAntennaFlagger.cc +++ b/steps/test/unit/tAntennaFlagger.cc @@ -48,15 +48,15 @@ class TestInput final : public InputStep { constexpr int n_receivers = 12 * 48; std::vector<int> ant1(n_baselines); std::vector<int> ant2(n_baselines); - int st1 = 0; - int st2 = 0; + int station1 = 0; + int station2 = 0; for (int i = 0; i < n_baselines; ++i) { - ant1[i] = st1; - ant2[i] = st2; - if (++st2 == n_receivers) { - st2 = 0; - if (++st1 == n_receivers) { - st1 = 0; + ant1[i] = station1; + ant2[i] = station2; + if (++station2 == n_receivers) { + station2 = 0; + if (++station1 == n_receivers) { + station1 = 0; } } } -- GitLab From 1af9bdbc1764c001e84ae58fe5476d32cc3797b9 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Tue, 23 Aug 2022 16:04:22 +0200 Subject: [PATCH 58/61] Add comment --- antennaflagger/Flagger.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/antennaflagger/Flagger.cc b/antennaflagger/Flagger.cc index 63ee20ccd..b6e50b90f 100644 --- a/antennaflagger/Flagger.cc +++ b/antennaflagger/Flagger.cc @@ -144,6 +144,7 @@ void FindBadAntennas(size_t n_stations, size_t n_receivers_per_station, std::vector<std::vector<int>> bad_stations(n_stations); // Check XY, YX and YY + // TODO: find out why XX is not checked and/or make this configurable for (size_t cor = 1; cor < n_correlations; ++cor) { // The check sum is the list of statistics for the current correlation // - in: T2 check_sum[n_correlations][n_receivers] -- GitLab From 4dd30e7e2f74ff9e9f4aec0844a9e12b094baab9 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Tue, 23 Aug 2022 16:06:02 +0200 Subject: [PATCH 59/61] Prevent 0 argument in std::log --- antennaflagger/Flagger.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/antennaflagger/Flagger.cc b/antennaflagger/Flagger.cc index b6e50b90f..6bc2f583e 100644 --- a/antennaflagger/Flagger.cc +++ b/antennaflagger/Flagger.cc @@ -158,8 +158,8 @@ void FindBadAntennas(size_t n_stations, size_t n_receivers_per_station, for (size_t i = 0; i < n_receivers_per_station; ++i) { const size_t station_index = station * n_receivers_per_station + i; station_indices[i] = station_index; - T real = std::log(check_sum[station_index].real()); - T imag = std::log(check_sum[station_index].imag()); + T real = std::log(static_cast<T>(1.0) + check_sum[station_index].real()); + T imag = std::log(static_cast<T>(1.0) + check_sum[station_index].imag()); check_sum_station[0 * n_receivers_per_station + i] = real; check_sum_station[1 * n_receivers_per_station + i] = imag; } -- GitLab From b3d84e3be20c64d357d63d9dca6f63014dd36c27 Mon Sep 17 00:00:00 2001 From: Bram Veenboer <bram.veenboer@gmail.com> Date: Tue, 23 Aug 2022 16:18:02 +0200 Subject: [PATCH 60/61] Rename some variables --- steps/AntennaFlagger.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index 5323c9cbf..04f76cb43 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -162,11 +162,11 @@ AntennaFlagger::AntennaFlagger(InputStep* input, name_ = prefix; if (do_detect_outliers_) { - std::string antennaSet(input->getInfo().antennaSet()); - std::vector<std::string> antennaNames(input->getInfo().antennaNames()); + std::string antenna_set(input->getInfo().antennaSet()); + std::vector<std::string> antenna_names(input->getInfo().antennaNames()); size_t n_receivers_per_station = 0; size_t n_stations = 0; - if (antennaNames[0].substr(0, 3) == "A12") { + if (antenna_names[0].substr(0, 3) == "A12") { // AARTFAAC-12: all antennas in a station are used as individual receivers n_receivers_per_station = 48; const size_t n_antennas = input->getInfo().nantenna(); -- GitLab From 755bc9fbc144d8093ab5ffbbe914daaa14f61791 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Offringa?= <offringa@astron.nl> Date: Tue, 23 Aug 2022 16:51:46 +0200 Subject: [PATCH 61/61] Formatting --- antennaflagger/Flagger.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/antennaflagger/Flagger.cc b/antennaflagger/Flagger.cc index 6bc2f583e..943a0d4e9 100644 --- a/antennaflagger/Flagger.cc +++ b/antennaflagger/Flagger.cc @@ -158,8 +158,10 @@ void FindBadAntennas(size_t n_stations, size_t n_receivers_per_station, for (size_t i = 0; i < n_receivers_per_station; ++i) { const size_t station_index = station * n_receivers_per_station + i; station_indices[i] = station_index; - T real = std::log(static_cast<T>(1.0) + check_sum[station_index].real()); - T imag = std::log(static_cast<T>(1.0) + check_sum[station_index].imag()); + T real = + std::log(static_cast<T>(1.0) + check_sum[station_index].real()); + T imag = + std::log(static_cast<T>(1.0) + check_sum[station_index].imag()); check_sum_station[0 * n_receivers_per_station + i] = real; check_sum_station[1 * n_receivers_per_station + i] = imag; } -- GitLab