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