From 943247addda12eb24ed417c21922038a543e0bcc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9=20Offringa?= <offringa@gmail.com>
Date: Wed, 16 Jun 2021 13:43:57 +0000
Subject: [PATCH] Faster sumtreshold algorithms

---
 .gitlab-ci.yml                                |   1 +
 CMakeLists.txt                                |  12 +-
 algorithms/antennaselector.cpp                |   6 +-
 algorithms/antennaselector.h                  |   3 +-
 algorithms/enums.h                            |  26 +-
 algorithms/sumthreshold.cpp                   |   5 +
 algorithms/sumthreshold.h                     |   1 +
 algorithms/sumthresholdmissing.cpp            | 122 ++++++++-
 algorithms/sumthresholdmissing.h              |  32 ++-
 algorithms/testsetgenerator.cpp               | 176 ++++++-------
 algorithms/testsetgenerator.h                 |  26 +-
 algorithms/thresholdconfig.cpp                |  10 +-
 .../controllers/histogrampagecontroller.cpp   |   2 +-
 doc/source/class_data.rst                     |   2 +
 lua/datawrapper.cpp                           |  24 +-
 rfigui/controllers/rfiguicontroller.cpp       |   7 +-
 rfigui/controllers/rfiguicontroller.h         |   3 +-
 rfigui/rfiguimenu.cpp                         |   2 +-
 rfigui/rfiguimenu.h                           |   4 +-
 rfigui/rfiguiwindow.cpp                       |  13 +-
 rfigui/rfiguiwindow.h                         |  29 ++-
 rfigui/simulatedialog.cpp                     | 102 ++++----
 rfigui/simulatedialog.h                       |   6 +-
 structures/timefrequencydata.cpp              |   6 +-
 test/algorithms/sumthresholdmissingtest.cpp   |  77 +++++-
 test/algorithms/sumthresholdtest.cpp          | 241 +++++++-----------
 test/algorithms/testtools.cpp                 | 177 +++++++++++++
 test/algorithms/testtools.h                   |  28 ++
 test/experiments/defaultstrategyspeedtest.cpp | 135 +++++++---
 test/interface/interfacetest.cpp              |   4 +-
 test/lua/defaultstrategytest.cpp              |  24 +-
 31 files changed, 899 insertions(+), 407 deletions(-)
 create mode 100644 test/algorithms/testtools.cpp
 create mode 100644 test/algorithms/testtools.h

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 9b706130..828e6700 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -24,6 +24,7 @@ before_script:
 
 aoflagger:
   script:
+    - scripts/run-clang-format.sh
     - mkdir build
     - cd build
     - cmake ../
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4dcf95ae..61ff80bf 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -92,7 +92,7 @@ if("${isSystemDir}" STREQUAL "-1")
    set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
 endif("${isSystemDir}" STREQUAL "-1")
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror=vla -DNDEBUG -funroll-loops -O3 -std=c++11")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror=vla -DNDEBUG -O3 -std=c++11")
 
 if(CMAKE_BUILD_TYPE STREQUAL "Debug")
   message(STATUS "Debug build selected: setting linking flag --no-undefined")
@@ -351,6 +351,15 @@ set_target_properties(aoflagger-bin PROPERTIES OUTPUT_NAME aoflagger)
 target_link_libraries(aoflagger-bin aoflagger-lib ${ALL_LIBRARIES})
 install (TARGETS aoflagger-bin DESTINATION bin)
 
+# A number of files perform the 'core' high-performance floating point
+# operations. In these files, NaNs are avoided and thus -ffast-math is
+# allowed. Note that visibilities can be NaN hence this can not be turned
+# on for all files.
+set_source_files_properties(
+	algorithms/sumthreshold.cpp
+	algorithms/sumthresholdmissing.cpp
+	PROPERTIES COMPILE_FLAGS "-ffast-math -funroll-loops")
+
 add_executable(runtests EXCLUDE_FROM_ALL
 	test/runtests.cpp
 	
@@ -377,6 +386,7 @@ add_executable(runtests EXCLUDE_FROM_ALL
 	test/algorithms/sumthresholdmissingtest.cpp
 	test/algorithms/sumthresholdtest.cpp
 	test/algorithms/thresholdtoolstest.cpp
+	test/algorithms/testtools.cpp
 		
 	test/structures/image2dtest.cpp
 	test/structures/timefrequencydatatest.cpp
diff --git a/algorithms/antennaselector.cpp b/algorithms/antennaselector.cpp
index 38b75c21..ff4f00ca 100644
--- a/algorithms/antennaselector.cpp
+++ b/algorithms/antennaselector.cpp
@@ -17,7 +17,8 @@ std::vector<size_t> AntennaSelector::Run(
   for (size_t p = 0; p != 4; ++p) {
     double meanStddev = 0.0;
     stddevs.clear();
-    for (const std::pair<const size_t, DefaultStatistics>& antenna : antStatistics) {
+    for (const std::pair<const size_t, DefaultStatistics>& antenna :
+         antStatistics) {
       double stddev = StatisticsDerivator::GetStatisticAmplitude(
           QualityTablesFormatter::StandardDeviationStatistic, antenna.second,
           p);
@@ -33,7 +34,8 @@ std::vector<size_t> AntennaSelector::Run(
 
     size_t index = 0;
     double limit = _threshold * stddevOfStddev;
-    for (const std::pair<const size_t, DefaultStatistics>& antenna : antStatistics) {
+    for (const std::pair<const size_t, DefaultStatistics>& antenna :
+         antStatistics) {
       if (std::fabs(stddevs[index] - meanStddev) > limit ||
           stddevs[index] == 0.0) {
         if (antenna.second.count[p] != 0) badAntennas.insert(antenna.first);
diff --git a/algorithms/antennaselector.h b/algorithms/antennaselector.h
index 22aa0d95..84ab0f24 100644
--- a/algorithms/antennaselector.h
+++ b/algorithms/antennaselector.h
@@ -13,7 +13,8 @@ class AntennaSelector {
  public:
   AntennaSelector() : _threshold(5.0) {}
 
-  std::vector<std::size_t> Run(const StatisticsCollection& statisticsCollection);
+  std::vector<std::size_t> Run(
+      const StatisticsCollection& statisticsCollection);
 
  private:
   void addStatistic(unsigned antIndex, const DefaultStatistics& stats,
diff --git a/algorithms/enums.h b/algorithms/enums.h
index aeec1e60..f6a2db3b 100644
--- a/algorithms/enums.h
+++ b/algorithms/enums.h
@@ -3,10 +3,10 @@
 
 namespace algorithms {
 
-enum class RFITestSet
-{
+enum class RFITestSet {
   Empty,
   SpectralLines,
+  GaussianSpectralLines,
   FullBandBursts,
   HalfBandBursts,
   VaryingBursts,
@@ -20,12 +20,13 @@ enum class RFITestSet
   PolarizedSpike
 };
 constexpr inline RFITestSet RFITestSetFirst() { return RFITestSet::Empty; }
-constexpr inline RFITestSet RFITestSetLast() { return RFITestSet::PolarizedSpike; }
+constexpr inline RFITestSet RFITestSetLast() {
+  return RFITestSet::PolarizedSpike;
+}
 RFITestSet inline operator++(RFITestSet& x) {
-    return x = (RFITestSet)(std::underlying_type<RFITestSet>::type(x) + 1); 
+  return x = (RFITestSet)(std::underlying_type<RFITestSet>::type(x) + 1);
 }
-enum class BackgroundTestSet
-{
+enum class BackgroundTestSet {
   Empty,
   LowFrequency,
   HighFrequency,
@@ -37,12 +38,17 @@ enum class BackgroundTestSet
   FaintVariableSidelobeSource,
   Checker
 };
-constexpr inline BackgroundTestSet BackgroundTestSetFirst() { return BackgroundTestSet::Empty; }
-constexpr inline BackgroundTestSet BackgroundTestSetLast() { return BackgroundTestSet::Checker; }
+constexpr inline BackgroundTestSet BackgroundTestSetFirst() {
+  return BackgroundTestSet::Empty;
+}
+constexpr inline BackgroundTestSet BackgroundTestSetLast() {
+  return BackgroundTestSet::Checker;
+}
 BackgroundTestSet inline operator++(BackgroundTestSet& x) {
-    return x = (BackgroundTestSet)(std::underlying_type<BackgroundTestSet>::type(x) + 1); 
+  return x = (BackgroundTestSet)(
+             std::underlying_type<BackgroundTestSet>::type(x) + 1);
 }
 
-}
+}  // namespace algorithms
 
 #endif
diff --git a/algorithms/sumthreshold.cpp b/algorithms/sumthreshold.cpp
index ef4557dc..7d5b2eec 100644
--- a/algorithms/sumthreshold.cpp
+++ b/algorithms/sumthreshold.cpp
@@ -16,6 +16,11 @@
 
 namespace algorithms {
 
+SumThreshold::VerticalScratch::VerticalScratch()
+    : lastFlaggedPos(nullptr, [](void*) {}),
+      sum(nullptr, [](void*) {}),
+      count(nullptr, [](void*) {}) {}
+
 SumThreshold::VerticalScratch::VerticalScratch(size_t width, size_t)
     : lastFlaggedPos(static_cast<int*>(boost::alignment::aligned_alloc(
                          64, sizeof(int) * width)),
diff --git a/algorithms/sumthreshold.h b/algorithms/sumthreshold.h
index 0265e702..2504ccee 100644
--- a/algorithms/sumthreshold.h
+++ b/algorithms/sumthreshold.h
@@ -12,6 +12,7 @@ namespace algorithms {
 class SumThreshold {
  public:
   struct VerticalScratch {
+    VerticalScratch();
     VerticalScratch(size_t width, size_t height);
     std::unique_ptr<int[], decltype(&free)> lastFlaggedPos;
     std::unique_ptr<num_t[], decltype(&free)> sum;
diff --git a/algorithms/sumthresholdmissing.cpp b/algorithms/sumthresholdmissing.cpp
index 13c42a83..b38813aa 100644
--- a/algorithms/sumthresholdmissing.cpp
+++ b/algorithms/sumthresholdmissing.cpp
@@ -2,6 +2,10 @@
 
 #include "../../structures/xyswappedmask2d.h"
 
+#include <vector>
+
+#include "sumthreshold.h"
+
 namespace algorithms {
 
 template <typename ImageLike, typename MaskLike, typename CMaskLike>
@@ -40,8 +44,8 @@ void SumThresholdMissing::horizontal(const ImageLike& input, MaskLike& mask,
           sum += input.Value(xRight, y);
           ++countAdded;
         }
-        // Check
-        if (countAdded > 0 && fabs(sum / countAdded) > threshold) {
+        // Threshold
+        if (countAdded > 0 && std::fabs(sum / countAdded) > threshold) {
           scratch.SetHorizontalValues(xLeft, y, true, xRight - xLeft + 1);
         }
         // subtract one sample at the left
@@ -69,9 +73,117 @@ template void SumThresholdMissing::horizontal(const Image2D& input,
                                               Mask2D& scratch, size_t length,
                                               num_t threshold);
 
-void SumThresholdMissing::Vertical(const Image2D& input, Mask2D& mask,
-                                   const Mask2D& missing, Mask2D& scratch,
-                                   size_t length, num_t threshold) {
+struct RowData {
+  num_t sum = 0.0f;
+  size_t yStart = 0;
+  size_t nNotFlagged = 0;
+  size_t nNotMissing = 0;
+};
+
+void SumThresholdMissing::InitializeVertical(VerticalCache& cache,
+                                             const Image2D& input,
+                                             const Mask2D& missing) {
+  cache.positions.assign(input.Width(), 0);
+  cache.validImage = Image2D::MakeSetImage(input.Width(), input.Height(), 0.0f);
+  cache.validMask = Mask2D::MakeSetMask<false>(input.Width(), input.Height());
+  for (size_t y = 0; y != input.Height(); ++y) {
+    for (size_t x = 0; x != input.Width(); ++x) {
+      if (!missing.Value(x, y)) {
+        size_t& pos = cache.positions[x];
+        cache.validImage.SetValue(x, pos, input.Value(x, y));
+        ++pos;
+      }
+    }
+  }
+  cache.scratch = SumThreshold::VerticalScratch(input.Width(), input.Height());
+}
+
+void SumThresholdMissing::VerticalStacked(VerticalCache& cache,
+                                          const Image2D& input, Mask2D& mask,
+                                          const Mask2D& missing,
+                                          Mask2D& scratch, size_t length,
+                                          num_t threshold) {
+  cache.positions.assign(cache.positions.size(), 0);
+  for (size_t y = 0; y != input.Height(); ++y) {
+    for (size_t x = 0; x != input.Width(); ++x) {
+      if (!missing.Value(x, y)) {
+        size_t& pos = cache.positions[x];
+        cache.validMask.SetValue(x, pos, mask.Value(x, y));
+        ++pos;
+      }
+    }
+  }
+
+  SumThreshold::VerticalLarge(&cache.validImage, &cache.validMask, &scratch,
+                              &cache.scratch, length, threshold);
+
+  cache.positions.assign(cache.positions.size(), 0);
+  for (size_t y = 0; y != input.Height(); ++y) {
+    for (size_t x = 0; x != input.Width(); ++x) {
+      size_t& pos = cache.positions[x];
+      if (!missing.Value(x, y)) {
+        if (cache.validMask.Value(x, pos)) mask.SetValue(x, y, true);
+        ++pos;
+      }
+    }
+  }
+}
+
+void SumThresholdMissing::VerticalConsecutive(const Image2D& input,
+                                              Mask2D& mask,
+                                              const Mask2D& missing,
+                                              Mask2D& scratch, size_t length,
+                                              num_t threshold) {
+  scratch = mask;
+  const size_t width = mask.Width(), height = mask.Height();
+  std::vector<RowData> rows(width);
+  if (length <= height) {
+    for (size_t y = 0; y != height; ++y) {
+      for (size_t x = 0; x != width; ++x) {
+        RowData& row = rows[x];
+
+        // Add sample if necessary
+        if (!missing.Value(x, y)) {
+          if (!mask.Value(x, y)) {
+            row.sum += input.Value(x, y);
+            ++row.nNotFlagged;
+          }
+          ++row.nNotMissing;
+          if (row.nNotMissing == length) {
+            // There is a full sequence after having added one more sample:
+            // perform threshold and subtract a sample that runs out
+            // of the window
+
+            // Threshold
+            if (row.nNotFlagged != 0 &&
+                std::fabs(row.sum / row.nNotFlagged) > threshold) {
+              scratch.SetVerticalValues(x, row.yStart, true,
+                                        y - row.yStart + 1);
+            }
+
+            // Subtract the oldest sample
+            do {
+              if (!missing.Value(x, row.yStart)) {
+                if (!mask.Value(x, row.yStart)) {
+                  row.sum -= input.Value(x, row.yStart);
+                  --row.nNotFlagged;
+                }
+                --row.nNotMissing;
+              }
+              ++row.yStart;
+            } while (row.nNotMissing == length);
+          }
+        }
+      }
+    }
+  }
+  mask = std::move(scratch);
+}
+
+void SumThresholdMissing::VerticalReference(const Image2D& input, Mask2D& mask,
+                                            const Mask2D& missing,
+                                            Mask2D& scratch, size_t length,
+                                            num_t threshold) {
   XYSwappedImage2D<const Image2D> swappedInput(input);
   XYSwappedMask2D<Mask2D> swappedMask(mask);
   XYSwappedMask2D<const Mask2D> swappedMissing(missing);
diff --git a/algorithms/sumthresholdmissing.h b/algorithms/sumthresholdmissing.h
index 43f40908..3ec3067f 100644
--- a/algorithms/sumthresholdmissing.h
+++ b/algorithms/sumthresholdmissing.h
@@ -4,19 +4,49 @@
 #include "../../structures/image2d.h"
 #include "../../structures/mask2d.h"
 
+#include "sumthreshold.h"
+
+#include <vector>
+
 namespace algorithms {
 
 class SumThresholdMissing {
  public:
+  struct VerticalCache {
+    std::vector<size_t> positions;
+    Image2D validImage;
+    Mask2D validMask;
+    SumThreshold::VerticalScratch scratch;
+  };
+
   static void Horizontal(const Image2D& input, Mask2D& mask,
                          const Mask2D& missing, Mask2D& scratch, size_t length,
                          num_t threshold) {
     horizontal(input, mask, missing, scratch, length, threshold);
   }
 
+  static void VerticalReference(const Image2D& input, Mask2D& mask,
+                                const Mask2D& missing, Mask2D& scratch,
+                                size_t length, num_t threshold);
+
+  static void InitializeVertical(VerticalCache& cache, const Image2D& input,
+                                 const Mask2D& missing);
+
   static void Vertical(const Image2D& input, Mask2D& mask,
                        const Mask2D& missing, Mask2D& scratch, size_t length,
-                       num_t threshold);
+                       num_t threshold) {
+    // VerticalReference(input, mask, missing, scratch, length, threshold);
+    // VerticalStacked(cache, input, mask, missing, scratch, length, threshold);
+    VerticalConsecutive(input, mask, missing, scratch, length, threshold);
+  }
+
+  static void VerticalConsecutive(const Image2D& input, Mask2D& mask,
+                                  const Mask2D& missing, Mask2D& scratch,
+                                  size_t length, num_t threshold);
+
+  static void VerticalStacked(VerticalCache& cache, const Image2D& input,
+                              Mask2D& mask, const Mask2D& missing,
+                              Mask2D& scratch, size_t length, num_t threshold);
 
  private:
   template <typename ImageLike, typename MaskLike, typename CMaskLike>
diff --git a/algorithms/testsetgenerator.cpp b/algorithms/testsetgenerator.cpp
index d6e781bb..8da1d447 100644
--- a/algorithms/testsetgenerator.cpp
+++ b/algorithms/testsetgenerator.cpp
@@ -21,22 +21,23 @@
 
 namespace algorithms {
 
-void TestSetGenerator::AddSpectralLine(Image2D& data, Mask2D& rfi, double lineStrength,
-                          size_t startChannel, size_t nChannels,
-                          double timeRatio,
-                          double timeOffsetRatio)
-{
+void TestSetGenerator::AddSpectralLine(Image2D& data, Mask2D& rfi,
+                                       double lineStrength, size_t startChannel,
+                                       size_t nChannels, double timeRatio,
+                                       double timeOffsetRatio,
+                                       enum BroadbandShape shape) {
   const size_t width = data.Width();
   const size_t tStart = size_t(timeOffsetRatio * width);
   const size_t tEnd = size_t((timeOffsetRatio + timeRatio) * width);
-  //const double tDuration = tEnd - tStart;
+  const double tDuration = tEnd - tStart;
   for (size_t t = tStart; t != tEnd; ++t) {
     // x will run from -1 to 1
-    // const double x = (double)((t - tStart) * 2) / tDuration - 1.0;
-    // double factor = shapeLevel(shape, x);
+    const double x = (double)((t - tStart) * 2) / tDuration - 1.0;
+    double factor = shapeLevel(shape, x);
     for (size_t ch = startChannel; ch < startChannel + nChannels; ++ch) {
-      data.AddValue(t, ch, lineStrength);
-      if (lineStrength > 0.0) rfi.SetValue(t, ch, true);
+      double value = lineStrength * factor;
+      data.AddValue(t, ch, value);
+      if (value != 0.0) rfi.SetValue(t, ch, true);
     }
   }
 }
@@ -128,8 +129,7 @@ Image2D TestSetGenerator::MakeGaussianData(unsigned width, unsigned height) {
 }
 
 std::string TestSetGenerator::GetDescription(BackgroundTestSet backgroundSet) {
-  switch(backgroundSet)
-  {
+  switch (backgroundSet) {
     case BackgroundTestSet::Empty:
       return "Empty";
     case BackgroundTestSet::LowFrequency:
@@ -160,6 +160,8 @@ std::string TestSetGenerator::GetDescription(RFITestSet rfiSet) {
       return "Empty";
     case RFITestSet::SpectralLines:
       return "Spectral lines";
+    case RFITestSet::GaussianSpectralLines:
+      return "Gaussian spectral lines";
     case RFITestSet::FullBandBursts:
       return "Full-band bursts / A";
     case RFITestSet::HalfBandBursts:
@@ -186,23 +188,24 @@ std::string TestSetGenerator::GetDescription(RFITestSet rfiSet) {
   return std::string();
 }
 
-TimeFrequencyData TestSetGenerator::MakeTestSet(RFITestSet rfiSet, BackgroundTestSet backgroundSet, size_t width, size_t height)
-{
+TimeFrequencyData TestSetGenerator::MakeTestSet(RFITestSet rfiSet,
+                                                BackgroundTestSet backgroundSet,
+                                                size_t width, size_t height) {
   TimeFrequencyData data;
-  if(rfiSet == algorithms::RFITestSet::PolarizedSpike)
-  {
+  if (rfiSet == algorithms::RFITestSet::PolarizedSpike) {
     std::vector<Image2DPtr> images(8);
     images[0] = Image2D::CreateZeroImagePtr(width, height);
     for (size_t p = 1; p != 8; ++p) {
       images[p] = images[0];
     }
     data = TimeFrequencyData::FromLinear(images[0], images[1], images[2],
-                                        images[3], images[4], images[5],
-                                        images[6], images[7]);
-  }
-  else {
-    Image2DPtr real = Image2D::MakePtr(TestSetGenerator::MakeNoise(width, height, 1.0));
-    Image2DPtr imag = Image2D::MakePtr(TestSetGenerator::MakeNoise(width, height, 1.0));
+                                         images[3], images[4], images[5],
+                                         images[6], images[7]);
+  } else {
+    Image2DPtr real =
+        Image2D::MakePtr(TestSetGenerator::MakeNoise(width, height, 1.0));
+    Image2DPtr imag =
+        Image2D::MakePtr(TestSetGenerator::MakeNoise(width, height, 1.0));
     data = TimeFrequencyData(aocommon::Polarization::StokesI, real, imag);
   }
   TestSetGenerator::MakeTestSet(rfiSet, data);
@@ -210,12 +213,13 @@ TimeFrequencyData TestSetGenerator::MakeTestSet(RFITestSet rfiSet, BackgroundTes
   return data;
 }
 
-void TestSetGenerator::MakeBackground(BackgroundTestSet testSet, TimeFrequencyData& data)
-{
-  if(testSet == BackgroundTestSet::StaticSidelobeSource || testSet == BackgroundTestSet::StrongVariableSidelobeSource || testSet == BackgroundTestSet::FaintVariableSidelobeSource)
-  {
+void TestSetGenerator::MakeBackground(BackgroundTestSet testSet,
+                                      TimeFrequencyData& data) {
+  if (testSet == BackgroundTestSet::StaticSidelobeSource ||
+      testSet == BackgroundTestSet::StrongVariableSidelobeSource ||
+      testSet == BackgroundTestSet::FaintVariableSidelobeSource) {
     SourceSet source;
-    switch(testSet) {
+    switch (testSet) {
       default:
       case BackgroundTestSet::StaticSidelobeSource:
         source = SourceSet::ConstantDistortion;
@@ -231,45 +235,44 @@ void TestSetGenerator::MakeBackground(BackgroundTestSet testSet, TimeFrequencyDa
     double bandwidth;
     bandwidth = (double)channelCount / 16.0 * 2500000.0;
     std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> pair =
-        DefaultModels::LoadSet(DefaultModels::NCPSet, source,
-                               0.0, channelCount, bandwidth);
+        DefaultModels::LoadSet(DefaultModels::NCPSet, source, 0.0, channelCount,
+                               bandwidth);
     data = pair.first;
-  }
-  else {
-    for(size_t i=0; i!=data.PolarizationCount(); ++i)
-    {
+  } else {
+    for (size_t i = 0; i != data.PolarizationCount(); ++i) {
       aocommon::PolarizationEnum pol = data.GetPolarization(i);
       TimeFrequencyData polData = data.MakeFromPolarizationIndex(i);
-      for(size_t imageIndex=0; imageIndex!=polData.ImageCount(); ++imageIndex)
-      {
-        const bool isImaginary = imageIndex==1;
-        if(isImaginary) {
+      for (size_t imageIndex = 0; imageIndex != polData.ImageCount();
+           ++imageIndex) {
+        const bool isImaginary = imageIndex == 1;
+        if (isImaginary) {
           Image2DPtr image = Image2D::MakePtr(*polData.GetImage(imageIndex));
           const size_t width = image->Width();
           const size_t height = image->Height();
-          switch(testSet) {
+          switch (testSet) {
             default:
             case BackgroundTestSet::Empty:
               break;
             case BackgroundTestSet::LowFrequency:
               for (unsigned y = 0; y < image->Height(); ++y) {
                 for (unsigned x = 0; x < image->Width(); ++x) {
-                  image->AddValue(x, y,
-                                sinn(num_t(x) * M_PIn * 5.0 / image->Width()) + 0.1);
+                  image->AddValue(
+                      x, y,
+                      sinn(num_t(x) * M_PIn * 5.0 / image->Width()) + 0.1);
                 }
               }
               break;
             case BackgroundTestSet::HighFrequency:
               for (unsigned y = 0; y < image->Height(); ++y) {
                 for (unsigned x = 0; x < image->Width(); ++x) {
-                  image->AddValue(
-                      x, y,
-                      sinn((long double)(x + y * 0.1) * M_PIn * 5.0L / image->Width() +
-                          0.1));
                   image->AddValue(x, y,
-                                sinn((long double)(x + pown(y, 1.1)) * M_PIn * 50.0L /
-                                          image->Width() +
-                                      0.1));
+                                  sinn((long double)(x + y * 0.1) * M_PIn *
+                                           5.0L / image->Width() +
+                                       0.1));
+                  image->AddValue(x, y,
+                                  sinn((long double)(x + pown(y, 1.1)) * M_PIn *
+                                           50.0L / image->Width() +
+                                       0.1));
                 }
               }
               for (unsigned y = 0; y < image->Height(); ++y) {
@@ -289,12 +292,12 @@ void TestSetGenerator::MakeBackground(BackgroundTestSet testSet, TimeFrequencyDa
               SubtractBackground(*image);
               break;
             case BackgroundTestSet::Checker:
-            for (size_t y = 0; y != height; ++y) {
-              for (size_t x = 0; x != width; ++x) {
-                image->SetValue(x, y, ((x + y) % 2 == 0) ? -1 : 2);
+              for (size_t y = 0; y != height; ++y) {
+                for (size_t x = 0; x != width; ++x) {
+                  image->SetValue(x, y, ((x + y) % 2 == 0) ? -1 : 2);
+                }
               }
-            }
-            break;
+              break;
           }
           polData.SetImage(imageIndex, image);
         }
@@ -304,63 +307,65 @@ void TestSetGenerator::MakeBackground(BackgroundTestSet testSet, TimeFrequencyDa
   }
 }
 
-void TestSetGenerator::MakeTestSet(RFITestSet testSet, TimeFrequencyData& data) {
-  for(size_t i=0; i!=data.PolarizationCount(); ++i)
-  {
+void TestSetGenerator::MakeTestSet(RFITestSet testSet,
+                                   TimeFrequencyData& data) {
+  for (size_t i = 0; i != data.PolarizationCount(); ++i) {
     aocommon::PolarizationEnum pol = data.GetPolarization(i);
     TimeFrequencyData polData = data.MakeFromPolarizationIndex(i);
     Mask2DPtr rfi(Mask2D::MakePtr(*polData.GetSingleMask()));
-    for(size_t imageIndex=0; imageIndex!=polData.ImageCount(); ++imageIndex)
-    {
-      const bool isImaginary = imageIndex==1;
+    for (size_t imageIndex = 0; imageIndex != polData.ImageCount();
+         ++imageIndex) {
+      const bool isImaginary = imageIndex == 1;
       Image2DPtr image = Image2D::MakePtr(*polData.GetImage(imageIndex));
       const size_t width = image->Width();
       const size_t height = image->Height();
-      if(pol == 2 && isImaginary && testSet == RFITestSet::PolarizedSpike) {
+      if (pol == 2 && isImaginary && testSet == RFITestSet::PolarizedSpike) {
         // Set one absurdly high value which all strategies should detect
         size_t rfiX = width / 2, rfiY = height / 2;
         image->SetValue(rfiX, rfiY, 1000.0);
-      }
-      else if(!isImaginary) {
+      } else if (!isImaginary) {
         switch (testSet) {
           case RFITestSet::Empty:
-          break;
+            break;
           case RFITestSet::SpectralLines:
-            AddSpectralLinesToTestSet(*image, *rfi, 1.0);
-          break;
+            AddSpectralLinesToTestSet(*image, *rfi, 1.0, UniformShape);
+            break;
+          case RFITestSet::GaussianSpectralLines:
+            AddSpectralLinesToTestSet(*image, *rfi, 1.0, GaussianShape);
+            break;
           case RFITestSet::FullBandBursts:
             AddBroadbandToTestSet(*image, *rfi, 1.0);
-          break;
+            break;
           case RFITestSet::HalfBandBursts:
             AddBroadbandToTestSet(*image, *rfi, 0.5);
-          break;
+            break;
           case RFITestSet::VaryingBursts:
             AddVarBroadbandToTestSet(*image, *rfi);
-          break;
+            break;
           case RFITestSet::GaussianBursts:
             AddBroadbandToTestSet(*image, *rfi, 1.0, 1.0, GaussianShape);
-          break;
+            break;
           case RFITestSet::SinusoidalBursts:
             AddBroadbandToTestSet(*image, *rfi, 1.0, 1.0, SinusoidalShape);
-          break;
+            break;
           case RFITestSet::SlewedGaussians:
             AddSlewedBroadbandToTestSet(*image, *rfi, 1.0);
-          break;
+            break;
           case RFITestSet::FluctuatingBursts:
             AddBurstBroadbandToTestSet(*image, *rfi);
-          break;
-          case  RFITestSet::StrongPowerLaw:
+            break;
+          case RFITestSet::StrongPowerLaw:
             *image += sampleRFIDistribution(width, height, 1.0);
             rfi->SetAll<true>();
-          break;
-          case  RFITestSet::MediumPowerLaw:
+            break;
+          case RFITestSet::MediumPowerLaw:
             *image += sampleRFIDistribution(width, height, 0.1);
             rfi->SetAll<true>();
-          break;
-          case  RFITestSet::WeakPowerLaw:
+            break;
+          case RFITestSet::WeakPowerLaw:
             *image += sampleRFIDistribution(width, height, 0.01);
             rfi->SetAll<true>();
-          break;
+            break;
           case RFITestSet::PolarizedSpike:
             break;
         }
@@ -373,17 +378,17 @@ void TestSetGenerator::MakeTestSet(RFITestSet testSet, TimeFrequencyData& data)
 }
 
 void TestSetGenerator::AddSpectralLinesToTestSet(Image2D& image, Mask2D& rfi,
-                                             double baseStrength) {
-  for(size_t i=0; i!=10; ++i) {
+                                                 double baseStrength,
+                                                 enum BroadbandShape shape) {
+  for (size_t i = 0; i != 10; ++i) {
     const size_t channel = ((i * 2 + 1) * image.Height()) / 20;
-    const double strength = baseStrength * (1.0 + (i*2.0/10.0));
-    AddSpectralLine(image, rfi, strength, channel, 1, 1.0, 0.0);
+    const double strength = baseStrength * (1.0 + (i * 2.0 / 10.0));
+    AddSpectralLine(image, rfi, strength, channel, 1, 1.0, 0.0, shape);
   }
 }
 
 void TestSetGenerator::AddBroadbandToTestSet(Image2D& image, Mask2D& rfi,
-                                             double length,
-                                             double strength,
+                                             double length, double strength,
                                              enum BroadbandShape shape) {
   size_t frequencyCount = image.Height();
   unsigned step = image.Width() / 11;
@@ -462,9 +467,8 @@ void TestSetGenerator::AddVarBroadbandToTestSet(Image2D& image, Mask2D& rfi) {
   AddBroadbandLine(image, rfi, 1.6, step * 10, 1, 0.444377, 0.240526);
 }
 
-void TestSetGenerator::SetModelData(Image2D& image,
-                                    unsigned sources, size_t width,
-                                    size_t height) {
+void TestSetGenerator::SetModelData(Image2D& image, unsigned sources,
+                                    size_t width, size_t height) {
   class Model model;
   if (sources >= 5) {
     model.AddSource(0.1, 0.1, 0.5);
diff --git a/algorithms/testsetgenerator.h b/algorithms/testsetgenerator.h
index 16d21796..5fc04831 100644
--- a/algorithms/testsetgenerator.h
+++ b/algorithms/testsetgenerator.h
@@ -65,9 +65,9 @@ class TestSetGenerator {
                      frequencyRatio, 0.5L - frequencyRatio / 2.0L);
   }
   static void AddSpectralLine(Image2D& data, Mask2D& rfi, double lineStrength,
-                               size_t startChannel, size_t nChannels,
-                               double timeRatio,
-                               double timeOffsetRatio);
+                              size_t startChannel, size_t nChannels,
+                              double timeRatio, double timeOffsetRatio,
+                              enum BroadbandShape shape);
   static void AddBroadbandLine(Image2D& data, Mask2D& rfi, double lineStrength,
                                size_t startTime, size_t duration,
                                double frequencyRatio,
@@ -89,11 +89,17 @@ class TestSetGenerator {
 
   static std::string GetDescription(BackgroundTestSet backgroundSet);
   static std::string GetDescription(RFITestSet rfiSet);
-  static TimeFrequencyData MakeTestSet(RFITestSet rfiSet, BackgroundTestSet backgroundSet, size_t width, size_t height);
-  static void MakeBackground(BackgroundTestSet backgroundSet, TimeFrequencyData& image);
+  static TimeFrequencyData MakeTestSet(RFITestSet rfiSet,
+                                       BackgroundTestSet backgroundSet,
+                                       size_t width, size_t height);
+  static void MakeBackground(BackgroundTestSet backgroundSet,
+                             TimeFrequencyData& image);
   static void MakeTestSet(RFITestSet testSet, TimeFrequencyData& data);
+
  private:
-  static void AddSpectralLinesToTestSet(Image2D& image, Mask2D& rfi, double strength = 1.0);
+  static void AddSpectralLinesToTestSet(Image2D& image, Mask2D& rfi,
+                                        double strength,
+                                        enum BroadbandShape shape);
   static void AddGaussianBroadbandToTestSet(Image2D& image, Mask2D& rfi) {
     AddBroadbandToTestSet(image, rfi, 1.0, 1.0, GaussianShape);
   }
@@ -107,15 +113,15 @@ class TestSetGenerator {
     AddSlewedBroadbandToTestSet(image, rfi, 1.0);
   }
 
-  static void AddBroadbandToTestSet(Image2D& image, Mask2D& rfi,
-                                    double length, double strength = 1.0,
+  static void AddBroadbandToTestSet(Image2D& image, Mask2D& rfi, double length,
+                                    double strength = 1.0,
                                     enum BroadbandShape shape = UniformShape);
   static void AddSlewedBroadbandToTestSet(
       Image2D& image, Mask2D& rfi, double length, double strength = 1.0,
       double slewrate = 0.02, enum BroadbandShape shape = GaussianShape);
   static void AddVarBroadbandToTestSet(Image2D& image, Mask2D& rfi);
-  static void SetModelData(Image2D& image, unsigned sources,
-                           size_t width, size_t height);
+  static void SetModelData(Image2D& image, unsigned sources, size_t width,
+                           size_t height);
   static void SubtractBackground(Image2D& image);
   static Image2D sampleRFIDistribution(unsigned width, unsigned height,
                                        double ig_over_rsq);
diff --git a/algorithms/thresholdconfig.cpp b/algorithms/thresholdconfig.cpp
index 0e09c271..4744d67e 100644
--- a/algorithms/thresholdconfig.cpp
+++ b/algorithms/thresholdconfig.cpp
@@ -145,7 +145,10 @@ void ThresholdConfig::ExecuteWithMissing(const Image2D* image, Mask2D* mask,
           ? _horizontalOperations.size()
           : _verticalOperations.size();
   SumThreshold::VerticalScratch normalScratch(mask->Width(), mask->Height());
-  for (unsigned i = 0; i < operationCount; ++i) {
+  SumThresholdMissing::VerticalCache vMissingCache;
+  if (_method == SumThreshold && missing)
+    SumThresholdMissing::InitializeVertical(vMissingCache, *image, *missing);
+  for (size_t i = 0; i < operationCount; ++i) {
     switch (_method) {
       case SumThreshold:
 
@@ -174,8 +177,9 @@ void ThresholdConfig::ExecuteWithMissing(const Image2D* image, Mask2D* mask,
                 _verticalOperations[i].length,
                 _verticalOperations[i].threshold * frequencyFactor);
           else
-            SumThresholdMissing::Vertical(
-                *image, *mask, *missing, scratch, _verticalOperations[i].length,
+            SumThresholdMissing::VerticalStacked(
+                vMissingCache, *image, *mask, *missing, scratch,
+                _verticalOperations[i].length,
                 _verticalOperations[i].threshold * frequencyFactor);
         }
         break;
diff --git a/aoqplot/controllers/histogrampagecontroller.cpp b/aoqplot/controllers/histogrampagecontroller.cpp
index eaf08296..0331099b 100644
--- a/aoqplot/controllers/histogrampagecontroller.cpp
+++ b/aoqplot/controllers/histogrampagecontroller.cpp
@@ -285,7 +285,7 @@ void HistogramPageController::addHistogramToPlot(
       logxEnd = log10(i.binEnd());
     } else {
       logxStart = log10(x);
-      logxEnd = 0.0; // unused, but to avoid warning
+      logxEnd = 0.0;  // unused, but to avoid warning
     }
     if (derivative) {
       const double cslope = histogram.NormalizedSlope(x / deltaS, x * deltaS);
diff --git a/doc/source/class_data.rst b/doc/source/class_data.rst
index 2ca12923..a377f614 100644
--- a/doc/source/class_data.rst
+++ b/doc/source/class_data.rst
@@ -35,9 +35,11 @@ Method summary
     - :meth:`~Data.copy`
     - :meth:`~Data.flag_nans`
     - :meth:`~Data.flag_zeros`
+    - :meth:`~Data.invert_mask`
     - :meth:`~Data.join_mask`
     - :meth:`~Data.set_mask`
     - :meth:`~Data.set_mask_for_channel_range`
+    - :meth:`~Data.set_masked_visibilities`
     - :meth:`~Data.set_polarization_data`
     - :meth:`~Data.set_visibilities`
     - :meth:`~Data.__sub`
diff --git a/lua/datawrapper.cpp b/lua/datawrapper.cpp
index 2aedd476..1e8468df 100644
--- a/lua/datawrapper.cpp
+++ b/lua/datawrapper.cpp
@@ -294,8 +294,7 @@ int Data::has_metadata(lua_State* L) {
 int Data::invert_mask(lua_State* L) {
   aoflagger_lua::Data* data = reinterpret_cast<aoflagger_lua::Data*>(
       luaL_checkudata(L, 1, "AOFlaggerData"));
-  for(size_t i=0; i!=data->TFData().MaskCount(); ++i)
-  {
+  for (size_t i = 0; i != data->TFData().MaskCount(); ++i) {
     Mask2DPtr mask = Mask2D::MakePtr(*data->TFData().GetMask(i));
     mask->Invert();
     data->TFData().SetMask(i, std::move(mask));
@@ -375,26 +374,27 @@ int Data::set_masked_visibilities(lua_State* L) {
       luaL_checkudata(L, 1, "AOFlaggerData"));
   aoflagger_lua::Data* rhs = reinterpret_cast<aoflagger_lua::Data*>(
       luaL_checkudata(L, 2, "AOFlaggerData"));
-  if (rhs->TFData().ImageCount() != lhs->TFData().ImageCount() || rhs->TFData().PolarizationCount() != lhs->TFData().PolarizationCount()) {
+  if (rhs->TFData().ImageCount() != lhs->TFData().ImageCount() ||
+      rhs->TFData().PolarizationCount() != lhs->TFData().PolarizationCount()) {
     std::string err =
-        "set_masked_visibilities() was executed with inconsistent data types: right "
+        "set_masked_visibilities() was executed with inconsistent data types: "
+        "right "
         "hand side had " +
         std::to_string(rhs->TFData().ImageCount()) + ", destination had " +
         std::to_string(lhs->TFData().ImageCount());
     return luaL_error(L, err.c_str());
   }
-  for (size_t p = 0; p != rhs->TFData().PolarizationCount(); ++p)
-  {
+  for (size_t p = 0; p != rhs->TFData().PolarizationCount(); ++p) {
     TimeFrequencyData lPolData = lhs->TFData().MakeFromPolarizationIndex(p);
-    const TimeFrequencyData rPolData = rhs->TFData().MakeFromPolarizationIndex(p);
+    const TimeFrequencyData rPolData =
+        rhs->TFData().MakeFromPolarizationIndex(p);
     const Mask2DCPtr mask(lPolData.GetSingleMask());
-    for (size_t i = 0; i != rhs->TFData().ImageCount(); ++i)
-    {
+    for (size_t i = 0; i != rhs->TFData().ImageCount(); ++i) {
       const Image2DCPtr source(rhs->TFData().GetImage(i));
       Image2DPtr dest(Image2D::MakePtr(*lhs->TFData().GetImage(i)));
-      for(size_t y=0; y!=source->Height(); ++y) {
-        for(size_t x=0; x!=source->Width(); ++x) {
-          if(mask->Value(x, y)) dest->SetValue(x, y, source->Value(x, y));
+      for (size_t y = 0; y != source->Height(); ++y) {
+        for (size_t x = 0; x != source->Width(); ++x) {
+          if (mask->Value(x, y)) dest->SetValue(x, y, source->Value(x, y));
         }
       }
       lPolData.SetImage(i, std::move(dest));
diff --git a/rfigui/controllers/rfiguicontroller.cpp b/rfigui/controllers/rfiguicontroller.cpp
index 9afce69d..0276a52c 100644
--- a/rfigui/controllers/rfiguicontroller.cpp
+++ b/rfigui/controllers/rfiguicontroller.cpp
@@ -320,7 +320,9 @@ std::vector<std::string> RFIGuiController::RecentFiles() const {
   return shownFiles;
 }
 
-void RFIGuiController::OpenTestSet(algorithms::RFITestSet rfiSet, algorithms::BackgroundTestSet backgroundSet) {
+void RFIGuiController::OpenTestSet(
+    algorithms::RFITestSet rfiSet,
+    algorithms::BackgroundTestSet backgroundSet) {
   size_t width = 1024, height = 512;
   TimeFrequencyMetaDataCPtr metaData;
   if (IsImageLoaded()) {
@@ -333,7 +335,8 @@ void RFIGuiController::OpenTestSet(algorithms::RFITestSet rfiSet, algorithms::Ba
   }
   CloseImageSet();
 
-  TimeFrequencyData data = TestSetGenerator::MakeTestSet(rfiSet, backgroundSet, width, height);
+  TimeFrequencyData data =
+      TestSetGenerator::MakeTestSet(rfiSet, backgroundSet, width, height);
   _tfController.SetNewData(data, metaData);
   const char* name = "Simulated test set";
   _tfController.Plot().SetTitleText(name);
diff --git a/rfigui/controllers/rfiguicontroller.h b/rfigui/controllers/rfiguicontroller.h
index e3e00506..4161a324 100644
--- a/rfigui/controllers/rfiguicontroller.h
+++ b/rfigui/controllers/rfiguicontroller.h
@@ -119,7 +119,8 @@ class RFIGuiController {
   void OpenMS(const std::vector<std::string>& filenames,
               const class MSOptions& options);
 
-  void OpenTestSet(algorithms::RFITestSet rfiSet, algorithms::BackgroundTestSet backgroundSet);
+  void OpenTestSet(algorithms::RFITestSet rfiSet,
+                   algorithms::BackgroundTestSet backgroundSet);
 
   std::vector<std::string> RecentFiles() const;
 
diff --git a/rfigui/rfiguimenu.cpp b/rfigui/rfiguimenu.cpp
index 28558a37..d5b9cf98 100644
--- a/rfigui/rfiguimenu.cpp
+++ b/rfigui/rfiguimenu.cpp
@@ -210,7 +210,7 @@ void RFIGuiMenu::makeBrowseMenu() {
 
 void RFIGuiMenu::makeSimulateMenu() {
   addItem(_menuSimulate, _miSimulate, OnSimulate, "Simulate...");
-  
+
   _miTestSetSubMenu.set_submenu(_menuTestSets);
   addItem(_menuSimulate, _miTestSetSubMenu, "Open _testset");
 
diff --git a/rfigui/rfiguimenu.h b/rfigui/rfiguimenu.h
index 17105954..ffbc1987 100644
--- a/rfigui/rfiguimenu.h
+++ b/rfigui/rfiguimenu.h
@@ -296,8 +296,8 @@ class RFIGuiMenu {
       _menuSimulate, _menuData, _menuSelectComplex, _menuSelectStokes,
       _menuSelectCircular, _menuSelectLinear, _menuSegmentation,
       _menuRecentFiles;
-  Gtk::MenuItem _miFile, _miView, _miStrategy, _miPlot, _miBrowse, _miSimulateMenu,
-      _miData;
+  Gtk::MenuItem _miFile, _miView, _miStrategy, _miPlot, _miBrowse,
+      _miSimulateMenu, _miData;
   ImgMenuItem _miRecentFiles;
 
   // File menu
diff --git a/rfigui/rfiguiwindow.cpp b/rfigui/rfiguiwindow.cpp
index 6f758891..eea5fafd 100644
--- a/rfigui/rfiguiwindow.cpp
+++ b/rfigui/rfiguiwindow.cpp
@@ -525,7 +525,8 @@ void RFIGuiWindow::onExecuteStrategyError(const std::string& error) {
   dialog.run();
 }
 
-void RFIGuiWindow::openTestSet(algorithms::RFITestSet rfiSet, algorithms::BackgroundTestSet backgroundSet) {
+void RFIGuiWindow::openTestSet(algorithms::RFITestSet rfiSet,
+                               algorithms::BackgroundTestSet backgroundSet) {
   _controller->OpenTestSet(rfiSet, backgroundSet);
 }
 
@@ -998,12 +999,12 @@ void RFIGuiWindow::showError(const std::string& description) {
   dialog.run();
 }
 
-void RFIGuiWindow::onSimulate()
-{
+void RFIGuiWindow::onSimulate() {
   SimulateDialog simDialog;
-  if(simDialog.run() == Gtk::RESPONSE_OK)
-  {
-    _controller->TFController().SetNewData(simDialog.Make(), TimeFrequencyMetaDataPtr(new TimeFrequencyMetaData()));
+  if (simDialog.run() == Gtk::RESPONSE_OK) {
+    _controller->TFController().SetNewData(
+        simDialog.Make(),
+        TimeFrequencyMetaDataPtr(new TimeFrequencyMetaData()));
     const char* name = "Simulated test set";
     _controller->TFController().Plot().SetTitleText(name);
     SetBaselineInfo(true, false, name, name);
diff --git a/rfigui/rfiguiwindow.h b/rfigui/rfiguiwindow.h
index ba62633f..509c79b1 100644
--- a/rfigui/rfiguiwindow.h
+++ b/rfigui/rfiguiwindow.h
@@ -157,11 +157,26 @@ class RFIGuiWindow : public Gtk::Window {
   void onOpenTestSetA() { openTestSet(algorithms::RFITestSet::FullBandBursts); }
   void onOpenTestSetB() { openTestSet(algorithms::RFITestSet::HalfBandBursts); }
   void onOpenTestSetC() { openTestSet(algorithms::RFITestSet::VaryingBursts); }
-  void onOpenTestSetD() { openTestSet(algorithms::RFITestSet::VaryingBursts, algorithms::BackgroundTestSet::ThreeSources); }
-  void onOpenTestSetE() { openTestSet(algorithms::RFITestSet::FullBandBursts, algorithms::BackgroundTestSet::FiveSources); }
-  void onOpenTestSetF() { openTestSet(algorithms::RFITestSet::VaryingBursts, algorithms::BackgroundTestSet::FiveSources); }
-  void onOpenTestSetG() { openTestSet(algorithms::RFITestSet::VaryingBursts, algorithms::BackgroundTestSet::FiveFilteredSources); }
-  void onOpenTestSetH() { openTestSet(algorithms::RFITestSet::VaryingBursts, algorithms::BackgroundTestSet::HighFrequency); }
+  void onOpenTestSetD() {
+    openTestSet(algorithms::RFITestSet::VaryingBursts,
+                algorithms::BackgroundTestSet::ThreeSources);
+  }
+  void onOpenTestSetE() {
+    openTestSet(algorithms::RFITestSet::FullBandBursts,
+                algorithms::BackgroundTestSet::FiveSources);
+  }
+  void onOpenTestSetF() {
+    openTestSet(algorithms::RFITestSet::VaryingBursts,
+                algorithms::BackgroundTestSet::FiveSources);
+  }
+  void onOpenTestSetG() {
+    openTestSet(algorithms::RFITestSet::VaryingBursts,
+                algorithms::BackgroundTestSet::FiveFilteredSources);
+  }
+  void onOpenTestSetH() {
+    openTestSet(algorithms::RFITestSet::VaryingBursts,
+                algorithms::BackgroundTestSet::HighFrequency);
+  }
   void onAddStaticFringe();
   void onAdd1SigmaFringe();
   void onSetToOne();
@@ -200,7 +215,9 @@ class RFIGuiWindow : public Gtk::Window {
   void onSimulate();
   void onHelpAbout();
 
-  void openTestSet(algorithms::RFITestSet rfiSet, algorithms::BackgroundTestSet backgroundSet = algorithms::BackgroundTestSet::Empty);
+  void openTestSet(algorithms::RFITestSet rfiSet,
+                   algorithms::BackgroundTestSet backgroundSet =
+                       algorithms::BackgroundTestSet::Empty);
 
   void onControllerStateChange();
 
diff --git a/rfigui/simulatedialog.cpp b/rfigui/simulatedialog.cpp
index 18787ec9..6bb62b38 100644
--- a/rfigui/simulatedialog.cpp
+++ b/rfigui/simulatedialog.cpp
@@ -6,16 +6,15 @@ using algorithms::BackgroundTestSet;
 using algorithms::RFITestSet;
 using algorithms::TestSetGenerator;
 
-SimulateDialog::SimulateDialog() :
-  _nTimesLabel("Number of timesteps:"),
-  _nChannelsLabel("Number of channels:"),
-  _bandwidthLabel("Bandwidth (MHz):"),
-  _polarizationsLabel("Polarizations:"),
-  _targetLabel("Target:"),
-  _noiseLabel("Noise:"),
-  _noiseLevelLabel("Noise level:"),
-  _rfiLabel("RFI:")
-{
+SimulateDialog::SimulateDialog()
+    : _nTimesLabel("Number of timesteps:"),
+      _nChannelsLabel("Number of channels:"),
+      _bandwidthLabel("Bandwidth (MHz):"),
+      _polarizationsLabel("Polarizations:"),
+      _targetLabel("Target:"),
+      _noiseLabel("Noise:"),
+      _noiseLevelLabel("Noise level:"),
+      _rfiLabel("RFI:") {
   _grid.attach(_nTimesLabel, 1, 0, 1, 1);
   _nTimesEntry.set_text("1000");
   _grid.attach(_nTimesEntry, 2, 0, 1, 1);
@@ -25,93 +24,100 @@ SimulateDialog::SimulateDialog() :
   _grid.attach(_bandwidthLabel, 1, 2, 1, 1);
   _bandwidthEntry.set_text("32");
   _grid.attach(_bandwidthEntry, 2, 2, 1, 1);
- 
+
   _grid.attach(_polarizationsLabel, 0, 3, 1, 1);
   _polarizationsSelection.append("Stokes I");
   _polarizationsSelection.append("XX, YY");
   _polarizationsSelection.append("XX, XY, YX, YY");
   _polarizationsSelection.set_active(0);
   _grid.attach(_polarizationsSelection, 1, 3, 2, 1);
-  
+
   _grid.attach(_targetLabel, 0, 4, 1, 1);
   _targetSelection.append("Current");
-  for(BackgroundTestSet bSet = algorithms::BackgroundTestSetFirst(); bSet != algorithms::BackgroundTestSetLast(); ++bSet) {
+  for (BackgroundTestSet bSet = algorithms::BackgroundTestSetFirst();
+       bSet != algorithms::BackgroundTestSetLast(); ++bSet) {
     _targetSelection.append(TestSetGenerator::GetDescription(bSet));
   }
-  _targetSelection.append(TestSetGenerator::GetDescription(algorithms::BackgroundTestSetLast()));
+  _targetSelection.append(
+      TestSetGenerator::GetDescription(algorithms::BackgroundTestSetLast()));
   _targetSelection.set_active(1);
   _grid.attach(_targetSelection, 1, 4, 2, 1);
-  
+
   _grid.attach(_noiseLabel, 0, 5, 1, 1);
   _noiseSelection.append("None");
   _noiseSelection.append("Complex Gaussian noise");
   _noiseSelection.append("Rayleigh noise");
   _noiseSelection.set_active(1);
   _grid.attach(_noiseSelection, 1, 5, 2, 1);
-  
+
   _grid.attach(_noiseLevelLabel, 1, 6, 1, 1);
   _noiseLevelEntry.set_text("1");
   _grid.attach(_noiseLevelEntry, 2, 6, 1, 1);
 
   _grid.attach(_rfiLabel, 0, 7, 1, 1);
-  for(RFITestSet rSet = algorithms::RFITestSetFirst(); rSet != algorithms::RFITestSetLast(); ++rSet) {
+  for (RFITestSet rSet = algorithms::RFITestSetFirst();
+       rSet != algorithms::RFITestSetLast(); ++rSet) {
     _rfiSelection.append(TestSetGenerator::GetDescription(rSet));
   }
-  _rfiSelection.append(TestSetGenerator::GetDescription(algorithms::RFITestSetLast()));
+  _rfiSelection.append(
+      TestSetGenerator::GetDescription(algorithms::RFITestSetLast()));
   _rfiSelection.set_active(0);
   _grid.attach(_rfiSelection, 1, 7, 2, 1);
- 
+
   get_content_area()->add(_grid);
   _grid.show_all();
-  
+
   add_button("Cancel", Gtk::RESPONSE_CANCEL);
   _simulateButton = add_button("Simulate", Gtk::RESPONSE_OK);
 }
 
-TimeFrequencyData SimulateDialog::Make() const
-{
+TimeFrequencyData SimulateDialog::Make() const {
   const size_t width = atoi(_nTimesEntry.get_text().c_str());
   const size_t height = atoi(_nChannelsEntry.get_text().c_str());
   const double noiseLevel = atof(_noiseLevelEntry.get_text().c_str());
 
   BackgroundTestSet backgroundSet = BackgroundTestSet::Empty;
   const bool replaceBackground = _targetSelection.get_active_row_number() != 0;
-  if(replaceBackground)
-    backgroundSet = BackgroundTestSet(_targetSelection.get_active_row_number()-1);
+  if (replaceBackground)
+    backgroundSet =
+        BackgroundTestSet(_targetSelection.get_active_row_number() - 1);
   const RFITestSet rfiSet = RFITestSet(_rfiSelection.get_active_row_number());
-  
+
   TimeFrequencyData data;
   const bool isComplex = _noiseSelection.get_active_row_number() != 2;
   const bool addNoise = _noiseSelection.get_active_row_number() != 0;
   std::vector<aocommon::PolarizationEnum> polarizations;
-  if(_polarizationsSelection.get_active_row_number() == 0)
-    polarizations = { aocommon::Polarization::StokesI };
-  else if(_polarizationsSelection.get_active_row_number() == 1)
-    polarizations = { aocommon::Polarization::XX, aocommon::Polarization::YY };
+  if (_polarizationsSelection.get_active_row_number() == 0)
+    polarizations = {aocommon::Polarization::StokesI};
+  else if (_polarizationsSelection.get_active_row_number() == 1)
+    polarizations = {aocommon::Polarization::XX, aocommon::Polarization::YY};
   else
-    polarizations = { aocommon::Polarization::XX, aocommon::Polarization::XY, aocommon::Polarization::YX, aocommon::Polarization::YY };
-  for(const aocommon::PolarizationEnum pol : polarizations)
-  {
-    if(isComplex) {
-      if(addNoise) {
-        Image2DPtr real = Image2D::MakePtr(TestSetGenerator::MakeNoise(width, height, noiseLevel));
-        Image2DPtr imag = Image2D::MakePtr(TestSetGenerator::MakeNoise(width, height, noiseLevel));
-        data = TimeFrequencyData::MakeFromPolarizationCombination(data, TimeFrequencyData(pol, real, imag));
-      }
-      else {
+    polarizations = {aocommon::Polarization::XX, aocommon::Polarization::XY,
+                     aocommon::Polarization::YX, aocommon::Polarization::YY};
+  for (const aocommon::PolarizationEnum pol : polarizations) {
+    if (isComplex) {
+      if (addNoise) {
+        Image2DPtr real = Image2D::MakePtr(
+            TestSetGenerator::MakeNoise(width, height, noiseLevel));
+        Image2DPtr imag = Image2D::MakePtr(
+            TestSetGenerator::MakeNoise(width, height, noiseLevel));
+        data = TimeFrequencyData::MakeFromPolarizationCombination(
+            data, TimeFrequencyData(pol, real, imag));
+      } else {
         Image2DPtr zero = Image2D::CreateZeroImagePtr(width, height);
-        data = TimeFrequencyData::MakeFromPolarizationCombination(data, TimeFrequencyData(pol, zero, zero));
+        data = TimeFrequencyData::MakeFromPolarizationCombination(
+            data, TimeFrequencyData(pol, zero, zero));
       }
-    }
-    else {
-      Image2DPtr amp = Image2D::MakePtr(TestSetGenerator::MakeRayleighData(width, height));
-      data = TimeFrequencyData::MakeFromPolarizationCombination(data, TimeFrequencyData(TimeFrequencyData::AmplitudePart, pol, amp));
+    } else {
+      Image2DPtr amp =
+          Image2D::MakePtr(TestSetGenerator::MakeRayleighData(width, height));
+      data = TimeFrequencyData::MakeFromPolarizationCombination(
+          data, TimeFrequencyData(TimeFrequencyData::AmplitudePart, pol, amp));
     }
   }
-  
-  if(replaceBackground)
-    TestSetGenerator::MakeBackground(backgroundSet, data);
+
+  if (replaceBackground) TestSetGenerator::MakeBackground(backgroundSet, data);
   TestSetGenerator::MakeTestSet(rfiSet, data);
-  
+
   return data;
 }
diff --git a/rfigui/simulatedialog.h b/rfigui/simulatedialog.h
index b612f064..d5c61d2f 100644
--- a/rfigui/simulatedialog.h
+++ b/rfigui/simulatedialog.h
@@ -18,7 +18,7 @@ class SimulateDialog : public Gtk::Dialog {
   ~SimulateDialog() {}
 
   TimeFrequencyData Make() const;
-  
+
  private:
   void onSimulateClicked();
   void onCloseClicked();
@@ -41,10 +41,8 @@ class SimulateDialog : public Gtk::Dialog {
   Gtk::Entry _noiseLevelEntry;
   Gtk::Label _rfiLabel;
   Gtk::ComboBoxText _rfiSelection;
-  
-  Gtk::Button* _simulateButton;
 
+  Gtk::Button* _simulateButton;
 };
 
 #endif  // IMAGEPROPERTIESWINDOW_H
-
diff --git a/structures/timefrequencydata.cpp b/structures/timefrequencydata.cpp
index cd29129a..5e40f6d3 100644
--- a/structures/timefrequencydata.cpp
+++ b/structures/timefrequencydata.cpp
@@ -156,10 +156,8 @@ TimeFrequencyData TimeFrequencyData::MakeFromPolarizationCombination(
 
 TimeFrequencyData TimeFrequencyData::MakeFromPolarizationCombination(
     const TimeFrequencyData &first, const TimeFrequencyData &second) {
-  if(first.IsEmpty())
-    return second;
-  if(second.IsEmpty())
-    return first;
+  if (first.IsEmpty()) return second;
+  if (second.IsEmpty()) return first;
   if (first.ComplexRepresentation() != second.ComplexRepresentation())
     throw std::runtime_error(
         "Trying to create polarization combination from data with different "
diff --git a/test/algorithms/sumthresholdmissingtest.cpp b/test/algorithms/sumthresholdmissingtest.cpp
index f6cc5687..ca28570b 100644
--- a/test/algorithms/sumthresholdmissingtest.cpp
+++ b/test/algorithms/sumthresholdmissingtest.cpp
@@ -1,3 +1,5 @@
+#include "testtools.h"
+
 #include "../../structures/image2d.h"
 #include "../../structures/mask2d.h"
 #include "../../structures/timefrequencydata.h"
@@ -6,12 +8,13 @@
 #include "../../algorithms/testsetgenerator.h"
 #include "../../algorithms/thresholdconfig.h"
 
-#include "../testingtools/maskasserter.h"
-
 #include <boost/test/unit_test.hpp>
 
 using algorithms::SumThresholdMissing;
 
+using test_tools::compareHorizontalAlgorithm;
+using test_tools::compareVerticalAlgorithm;
+
 BOOST_AUTO_TEST_SUITE(masked_sumthreshold,
                       *boost::unit_test::label("algorithms"))
 
@@ -38,4 +41,74 @@ BOOST_AUTO_TEST_CASE(horizontal) {
   }
 }
 
+BOOST_AUTO_TEST_CASE(vertical) {
+  const unsigned width = 8, height = 8;
+  Mask2D mask = Mask2D::MakeSetMask<false>(width, height);
+  Mask2D missing = Mask2D::MakeSetMask<false>(width, height);
+  Mask2D scratch = Mask2D::MakeUnsetMask(width, height);
+  Image2D image = Image2D::MakeSetImage(width, height, 0.0);
+
+  for (size_t x = 0; x != width; ++x) {
+    image.SetValue(x, 3, 1.0);
+    image.SetValue(x, 4, 1.0);
+  }
+
+  SumThresholdMissing::Horizontal(image, mask, missing, scratch, 2, 0.8);
+
+  for (size_t x = 0; x != width; ++x) {
+    BOOST_CHECK(mask.Value(x, 3));
+    BOOST_CHECK(mask.Value(x, 4));
+    BOOST_CHECK(!mask.Value(x, 0));
+    BOOST_CHECK(!mask.Value(x, 2));
+    BOOST_CHECK(!mask.Value(x, 5));
+  }
+}
+
+BOOST_AUTO_TEST_CASE(compare_vertical_masked_reference) {
+  compareVerticalAlgorithm(
+      [](const Image2D* input, Mask2D* mask, Mask2D* scratch, size_t length,
+         num_t threshold) {
+        Image2D mInput;
+        Mask2D mMask;
+        Mask2D missing;
+        test_tools::introduceGap(*input, *mask, mInput, mMask, missing);
+        SumThresholdMissing::VerticalReference(mInput, mMask, missing, *scratch,
+                                               length, threshold);
+        test_tools::removeGap(*mask, mMask);
+      },
+      {});
+}
+
+BOOST_AUTO_TEST_CASE(compare_vertical_masked_consecutive) {
+  compareVerticalAlgorithm(
+      [](const Image2D* input, Mask2D* mask, Mask2D* scratch, size_t length,
+         num_t threshold) {
+        Image2D mInput;
+        Mask2D mMask;
+        Mask2D missing;
+        test_tools::introduceGap(*input, *mask, mInput, mMask, missing);
+        SumThresholdMissing::VerticalConsecutive(mInput, mMask, missing,
+                                                 *scratch, length, threshold);
+        test_tools::removeGap(*mask, mMask);
+      },
+      {});
+}
+
+BOOST_AUTO_TEST_CASE(compare_vertical_masked_stacked) {
+  compareVerticalAlgorithm(
+      [](const Image2D* input, Mask2D* mask, Mask2D* scratch, size_t length,
+         num_t threshold) {
+        Image2D mInput;
+        Mask2D mMask;
+        Mask2D missing;
+        test_tools::introduceGap(*input, *mask, mInput, mMask, missing);
+        SumThresholdMissing::VerticalCache cache;
+        SumThresholdMissing::InitializeVertical(cache, mInput, missing);
+        SumThresholdMissing::VerticalStacked(cache, mInput, mMask, missing,
+                                             *scratch, length, threshold);
+        test_tools::removeGap(*mask, mMask);
+      },
+      {});
+}
+
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/test/algorithms/sumthresholdtest.cpp b/test/algorithms/sumthresholdtest.cpp
index fcd081b7..d24b317e 100644
--- a/test/algorithms/sumthresholdtest.cpp
+++ b/test/algorithms/sumthresholdtest.cpp
@@ -1,3 +1,5 @@
+#include "testtools.h"
+
 #include "../../structures/image2d.h"
 #include "../../structures/mask2d.h"
 #include "../../structures/timefrequencydata.h"
@@ -5,93 +7,84 @@
 #include "../../algorithms/testsetgenerator.h"
 #include "../../algorithms/thresholdconfig.h"
 #include "../../algorithms/sumthreshold.h"
-
-#include "../testingtools/maskasserter.h"
+#include "../../algorithms/sumthresholdmissing.h"
 
 #include <boost/test/unit_test.hpp>
 
-using aocommon::Polarization;
-
-using algorithms::BackgroundTestSet;
-using algorithms::RFITestSet;
 using algorithms::SumThreshold;
-using algorithms::TestSetGenerator;
+using algorithms::SumThresholdMissing;
 using algorithms::ThresholdConfig;
 
-BOOST_AUTO_TEST_SUITE(sumthreshold, *boost::unit_test::label("algorithms"))
-
-#ifdef __SSE__
-BOOST_AUTO_TEST_CASE(vertical_sumthreshold_SSE) {
-  const unsigned width = 2048, height = 256;
-  Mask2D mask2 = Mask2D::MakeUnsetMask(width, height);
-  Mask2D scratch = Mask2D::MakeUnsetMask(width, height);
-  TimeFrequencyData data =
-    TestSetGenerator::MakeTestSet(RFITestSet::GaussianBursts, BackgroundTestSet::Empty, width, height);
-  Mask2D mask1(*data.GetSingleMask());
-  Image2DCPtr image = data.GetSingleImage();
-
-  ThresholdConfig config;
-  config.InitializeLengthsDefault(9);
-  num_t mode = image->GetMode();
-  config.InitializeThresholdsFromFirstThreshold(6.0 * mode,
-                                                ThresholdConfig::Rayleigh);
-  for (unsigned i = 0; i < 9; ++i) {
-    mask1.SetAll<false>();
-    mask2.SetAll<false>();
+using test_tools::compareHorizontalAlgorithm;
+using test_tools::compareVerticalAlgorithm;
 
-    const unsigned length = config.GetHorizontalLength(i);
-    const double threshold = config.GetHorizontalThreshold(i);
-
-    SumThreshold::VerticalLargeReference(image.get(), &mask1, &scratch, length,
-                                         threshold);
-    SumThreshold::VerticalLargeSSE(image.get(), &mask2, &scratch, length,
-                                   threshold);
+BOOST_AUTO_TEST_SUITE(sumthreshold, *boost::unit_test::label("algorithms"))
 
-    if (length != 32) {
-      BOOST_CHECK(mask1 == mask2);
-    }
-  }
+BOOST_AUTO_TEST_CASE(masked_vertical_reference) {
+  compareVerticalAlgorithm(
+      [](const Image2D* input, Mask2D* mask, Mask2D* scratch, size_t length,
+         num_t threshold) {
+        Mask2D missing =
+            Mask2D::MakeSetMask<false>(input->Width(), input->Height());
+        SumThresholdMissing::VerticalReference(*input, *mask, missing, *scratch,
+                                               length, threshold);
+      },
+      {});
 }
 
-BOOST_AUTO_TEST_CASE(horizontal_sumthreshold_SSE) {
-  const unsigned width = 2048, height = 256;
-  Mask2D mask1 = Mask2D::MakeUnsetMask(width, height),
-         mask2 = Mask2D::MakeUnsetMask(width, height),
-         scratch = Mask2D::MakeUnsetMask(width, height);
-  TimeFrequencyData data =
-    TestSetGenerator::MakeTestSet(RFITestSet::GaussianBursts, BackgroundTestSet::Empty, width, height);
-  Image2DPtr real(Image2D::MakePtr(*data.GetImage(0)));
-  Image2DPtr imag(Image2D::MakePtr(*data.GetImage(1)));
-
-  mask1.SwapXY();
-  mask2.SwapXY();
-  real->SwapXY();
-  imag->SwapXY();
-
-  Image2DCPtr image = TimeFrequencyData(Polarization::XX, real, imag).GetSingleImage();
+BOOST_AUTO_TEST_CASE(masked_vertical_consecutive) {
+  compareVerticalAlgorithm(
+      [](const Image2D* input, Mask2D* mask, Mask2D* scratch, size_t length,
+         num_t threshold) {
+        Mask2D missing =
+            Mask2D::MakeSetMask<false>(input->Width(), input->Height());
+        SumThresholdMissing::VerticalConsecutive(*input, *mask, missing,
+                                                 *scratch, length, threshold);
+      },
+      {});
+}
 
-  ThresholdConfig config;
-  config.InitializeLengthsDefault(9);
-  num_t mode = image->GetMode();
-  config.InitializeThresholdsFromFirstThreshold(6.0 * mode,
-                                                ThresholdConfig::Rayleigh);
-  for (unsigned i = 0; i < 9; ++i) {
-    mask1.SetAll<false>();
-    mask2.SetAll<false>();
+BOOST_AUTO_TEST_CASE(masked_vertical_stacked) {
+  compareVerticalAlgorithm(
+      [](const Image2D* input, Mask2D* mask, Mask2D* scratch, size_t length,
+         num_t threshold) {
+        Mask2D missing =
+            Mask2D::MakeSetMask<false>(input->Width(), input->Height());
+        SumThresholdMissing::VerticalCache cache;
+        SumThresholdMissing::InitializeVertical(cache, *input, missing);
+        SumThresholdMissing::VerticalStacked(cache, *input, *mask, missing,
+                                             *scratch, length, threshold);
+      },
+      {});
+}
 
-    const unsigned length = config.GetHorizontalLength(i);
-    const double threshold = config.GetHorizontalThreshold(i);
+BOOST_AUTO_TEST_CASE(masked_horizontal_reference) {
+  compareHorizontalAlgorithm([](const Image2D* input, Mask2D* mask,
+                                Mask2D* scratch, size_t length,
+                                num_t threshold) {
+    Mask2D missing =
+        Mask2D::MakeSetMask<false>(input->Width(), input->Height());
+    SumThresholdMissing::Horizontal(*input, *mask, missing, *scratch, length,
+                                    threshold);
+  });
+}
 
-    SumThreshold::HorizontalLargeReference(image.get(), &mask1, &scratch,
-                                           length, threshold);
-    SumThreshold::HorizontalLargeSSE(image.get(), &mask2, &scratch, length,
-                                     threshold);
+#ifdef __SSE__
+BOOST_AUTO_TEST_CASE(vertical_SSE) {
+  compareVerticalAlgorithm(
+      [](const Image2D* input, Mask2D* mask, Mask2D* scratch, size_t length,
+         num_t threshold) {
+        SumThreshold::VerticalLargeSSE(input, mask, scratch, length, threshold);
+      },
+      {32});
+}
 
-    std::stringstream s;
-    s << "Equal SSE and reference masks produced by SumThreshold length "
-      << length << ", threshold " << threshold;
-    MaskAsserter::AssertEqualMasks(mask2, mask1, s.str());
-  }
+BOOST_AUTO_TEST_CASE(horizontal_SSE) {
+  compareHorizontalAlgorithm([](const Image2D* input, Mask2D* mask,
+                                Mask2D* scratch, size_t length,
+                                num_t threshold) {
+    SumThreshold::HorizontalLargeSSE(input, mask, scratch, length, threshold);
+  });
 }
 
 BOOST_AUTO_TEST_CASE(stability_SSE) {
@@ -125,94 +118,32 @@ BOOST_AUTO_TEST_CASE(stability_SSE) {
 
 #ifdef __AVX2__
 BOOST_AUTO_TEST_CASE(horizontal_sumthreshold_AVX_dumas) {
-  const unsigned width = 2048, height = 256;
-  Mask2D mask1 = Mask2D::MakeUnsetMask(width, height),
-         mask2 = Mask2D::MakeUnsetMask(width, height),
-         scratch = Mask2D::MakeUnsetMask(width, height);
-  TimeFrequencyData data =
-    TestSetGenerator::MakeTestSet(RFITestSet::GaussianBursts, BackgroundTestSet::Empty, width, height);
-  Image2DPtr real(Image2D::MakePtr(*data.GetImage(0)));
-  Image2DPtr imag(Image2D::MakePtr(*data.GetImage(1)));
-
-  mask1.SwapXY();
-  mask2.SwapXY();
-  real->SwapXY();
-  imag->SwapXY();
-
-  Image2DCPtr image = TimeFrequencyData(Polarization::XX, real, imag).GetSingleImage();
-
-  ThresholdConfig config;
-  config.InitializeLengthsDefault(9);
-  num_t mode = image->GetMode();
-  config.InitializeThresholdsFromFirstThreshold(6.0 * mode,
-                                                ThresholdConfig::Rayleigh);
-  for (unsigned i = 0; i < 9; ++i) {
-    mask1.SetAll<false>();
-    mask2.SetAll<false>();
-
-    const unsigned length = config.GetHorizontalLength(i);
-    const double threshold = config.GetHorizontalThreshold(i);
-
-    SumThreshold::HorizontalLargeReference(image.get(), &mask1, &scratch,
-                                           length, threshold);
-    SumThreshold::HorizontalAVXDumas(image.get(), &mask2, length, threshold);
-
-    std::stringstream s;
-    s << "Equal AVX Dumas and reference masks produced by SumThreshold length "
-      << length << ", threshold " << threshold;
-    MaskAsserter::AssertEqualMasks(mask2, mask1, s.str());
-  }
-}
-
-static void VerticalSumThresholdAVX(bool dumas) {
-  const unsigned width = 2048, height = 256;
-  Mask2D mask1 = Mask2D::MakeUnsetMask(width, height),
-         mask2 = Mask2D::MakeUnsetMask(width, height),
-         scratch = Mask2D::MakeUnsetMask(width, height);
-  SumThreshold::VerticalScratch vScratch(width, height);
-  TimeFrequencyData data =
-    TestSetGenerator::MakeTestSet(RFITestSet::GaussianBursts, BackgroundTestSet::Empty, width, height);
-  Image2DCPtr image = data.GetSingleImage();
-
-  ThresholdConfig config;
-  config.InitializeLengthsDefault(9);
-  num_t mode = image->GetMode();
-  config.InitializeThresholdsFromFirstThreshold(6.0 * mode,
-                                                ThresholdConfig::Rayleigh);
-  for (unsigned i = 0; i < 9; ++i) {
-    mask1.SetAll<false>();
-    mask2.SetAll<false>();
-
-    const unsigned length = config.GetHorizontalLength(i);
-    const double threshold = config.GetHorizontalThreshold(i);
-
-    SumThreshold::VerticalLargeReference(image.get(), &mask1, &scratch, length,
-                                         threshold);
-    // std::cout << mask1.GetCount<true>() << '\n';
-    if (dumas)
-      SumThreshold::VerticalAVXDumas(image.get(), &mask2, &vScratch, length,
-                                     threshold);
-    else
-      SumThreshold::VerticalLargeAVX(image.get(), &mask2, &scratch, length,
-                                     threshold);
-
-    if (length != 32 || dumas) {  // I think there was a numerical issue on some
-                                  // platforms with 32, so just skip.
-      std::stringstream s;
-      s << "Equal AVX Dumas and reference masks produced by SumThreshold "
-           "length "
-        << length;
-      MaskAsserter::AssertEqualMasks(mask2, mask1, s.str());
-    }
-  }
+  compareHorizontalAlgorithm([](const Image2D* input, Mask2D* mask, Mask2D*,
+                                size_t length, num_t threshold) {
+    SumThreshold::HorizontalAVXDumas(input, mask, length, threshold);
+  });
 }
 
 BOOST_AUTO_TEST_CASE(vertical_sumthreshold_AVX) {
-  VerticalSumThresholdAVX(false);
+  compareVerticalAlgorithm(
+      [](const Image2D* input, Mask2D* mask, Mask2D*, size_t length,
+         num_t threshold) {
+        Mask2D scratch(*mask);
+        SumThreshold::VerticalLargeAVX(input, mask, &scratch, length,
+                                       threshold);
+      },
+      {32});
 }
 
 BOOST_AUTO_TEST_CASE(vertical_sumthreshold_AVX_dumas) {
-  VerticalSumThresholdAVX(true);
+  compareVerticalAlgorithm(
+      [](const Image2D* input, Mask2D* mask, Mask2D*, size_t length,
+         num_t threshold) {
+        SumThreshold::VerticalScratch vScratch(input->Width(), input->Height());
+        SumThreshold::VerticalAVXDumas(input, mask, &vScratch, length,
+                                       threshold);
+      },
+      {32});
 }
 
 static void SimpleVerticalSumThresholdAVX(bool dumas) {
diff --git a/test/algorithms/testtools.cpp b/test/algorithms/testtools.cpp
new file mode 100644
index 00000000..f2a7a17b
--- /dev/null
+++ b/test/algorithms/testtools.cpp
@@ -0,0 +1,177 @@
+#include "testtools.h"
+
+#include "../../structures/image2d.h"
+#include "../../structures/mask2d.h"
+#include "../../structures/timefrequencydata.h"
+
+#include "../../algorithms/testsetgenerator.h"
+#include "../../algorithms/thresholdconfig.h"
+#include "../../algorithms/sumthreshold.h"
+#include "../../algorithms/sumthresholdmissing.h"
+
+#include <boost/test/unit_test.hpp>
+
+#include <functional>
+#include <set>
+
+using aocommon::Polarization;
+
+using algorithms::BackgroundTestSet;
+using algorithms::RFITestSet;
+using algorithms::SumThreshold;
+using algorithms::TestSetGenerator;
+using algorithms::ThresholdConfig;
+
+namespace test_tools {
+
+namespace {
+constexpr size_t kFeatureSize = 16;
+}
+
+void compareVerticalAlgorithm(
+    std::function<void(const Image2D*, Mask2D*, Mask2D*, size_t, num_t)>
+        executeAlgorithm,
+    const std::set<size_t>& skips) {
+  const unsigned width = 512, height = 256;
+  Mask2D mask = Mask2D::MakeUnsetMask(width, height);
+  Mask2D scratch = Mask2D::MakeUnsetMask(width, height);
+  TimeFrequencyData data = TestSetGenerator::MakeTestSet(
+      RFITestSet::GaussianBursts, BackgroundTestSet::Empty, width, height);
+  Mask2D referenceMask(*data.GetSingleMask());
+  Image2DCPtr image = data.GetSingleImage();
+
+  ThresholdConfig config;
+  config.InitializeLengthsDefault(9);
+  num_t mode = image->GetMode();
+  config.InitializeThresholdsFromFirstThreshold(6.0 * mode,
+                                                ThresholdConfig::Rayleigh);
+  for (unsigned i = 0; i < 9; ++i) {
+    referenceMask.SetAll<false>();
+    mask.SetAll<false>();
+
+    const unsigned length = config.GetHorizontalLength(i);
+    const double threshold = config.GetHorizontalThreshold(i);
+
+    SumThreshold::VerticalLargeReference(image.get(), &referenceMask, &scratch,
+                                         length, threshold);
+    executeAlgorithm(image.get(), &mask, &scratch, length, threshold);
+
+    if (skips.find(length) == skips.end()) {
+      BOOST_CHECK_EQUAL(referenceMask.GetCount<true>(), mask.GetCount<true>());
+
+      std::ostringstream str;
+      str << "referenceMask == mask for length " << length;
+      BOOST_CHECK_MESSAGE(referenceMask == mask, str.str());
+    }
+  }
+}
+
+void compareHorizontalAlgorithm(
+    std::function<void(const Image2D*, Mask2D*, Mask2D*, size_t, num_t)>
+        executeAlgorithm) {
+  const size_t width = 512, height = 256;
+  Mask2D referenceMask = Mask2D::MakeUnsetMask(width, height),
+         mask = Mask2D::MakeUnsetMask(width, height),
+         scratch = Mask2D::MakeUnsetMask(width, height);
+  TimeFrequencyData data = TestSetGenerator::MakeTestSet(
+      RFITestSet::GaussianBursts, BackgroundTestSet::Empty, width, height);
+  Image2DPtr real(Image2D::MakePtr(*data.GetImage(0)));
+  Image2DPtr imag(Image2D::MakePtr(*data.GetImage(1)));
+
+  referenceMask.SwapXY();
+  mask.SwapXY();
+  real->SwapXY();
+  imag->SwapXY();
+
+  Image2DCPtr image =
+      TimeFrequencyData(Polarization::XX, real, imag).GetSingleImage();
+
+  ThresholdConfig config;
+  config.InitializeLengthsDefault(9);
+  num_t mode = image->GetMode();
+  config.InitializeThresholdsFromFirstThreshold(6.0 * mode,
+                                                ThresholdConfig::Rayleigh);
+  for (size_t i = 0; i < 9; ++i) {
+    referenceMask.SetAll<false>();
+    mask.SetAll<false>();
+
+    const size_t length = config.GetHorizontalLength(i);
+    const double threshold = config.GetHorizontalThreshold(i);
+
+    SumThreshold::HorizontalLargeReference(image.get(), &referenceMask,
+                                           &scratch, length, threshold);
+    executeAlgorithm(image.get(), &mask, &scratch, length, threshold);
+
+    std::ostringstream str;
+    str << "referenceMask == mask for length " << length;
+    BOOST_CHECK_MESSAGE(mask == referenceMask, str.str());
+  }
+}
+
+void introduceGap(const Image2D& input, const Mask2D& mask, Image2D& mInput,
+                  Mask2D& mMask, Mask2D& missing) {
+  missing = Mask2D::MakeSetMask<false>(input.Width() + kFeatureSize,
+                                       input.Height() + kFeatureSize);
+  mMask = Mask2D::MakeUnsetMask(mask.Width() + kFeatureSize,
+                                mask.Height() + kFeatureSize);
+  mInput = Image2D::MakeUnsetImage(input.Width() + kFeatureSize,
+                                   input.Height() + kFeatureSize);
+  const size_t x1 = input.Width() / 2;
+  const size_t x2 = x1 + kFeatureSize;
+  const size_t y1 = input.Height() / 2;
+  const size_t y2 = y1 + kFeatureSize;
+  for (size_t y = 0; y != y2; ++y) {
+    const num_t* dataRowIn = input.ValuePtr(0, y);
+    num_t* dataRowOut = mInput.ValuePtr(0, y);
+    std::copy_n(dataRowIn, x1, dataRowOut);
+    std::copy_n(dataRowIn, kFeatureSize, dataRowOut + x1);
+    std::copy_n(dataRowIn + x1, input.Width() - x1, dataRowOut + x2);
+    const bool* maskRowIn = mask.ValuePtr(0, y);
+    bool* maskRowOut = mMask.ValuePtr(0, y);
+    std::copy_n(maskRowIn, x1, maskRowOut);
+    std::copy_n(maskRowIn, kFeatureSize, maskRowOut + x1);
+    std::copy_n(maskRowIn + x1, input.Width() - x1, maskRowOut + x2);
+    std::fill_n(missing.ValuePtr(x1, y), kFeatureSize, true);
+  }
+  for (size_t y = y1; y != y2; ++y) {
+    const num_t* rowIn = input.ValuePtr(0, y);
+    num_t* rowOut = mInput.ValuePtr(0, y);
+    std::copy_n(rowIn, x1, rowOut);
+    std::copy_n(rowIn, kFeatureSize, rowOut + x1);
+    std::copy_n(rowIn + x1, input.Width() - x1, rowOut + x2);
+    std::fill_n(missing.ValuePtr(0, y), missing.Width(), true);
+  }
+  for (size_t y = y1; y != input.Height(); ++y) {
+    const num_t* dataRowIn = input.ValuePtr(0, y);
+    num_t* dataRowOut = mInput.ValuePtr(0, y + kFeatureSize);
+    std::copy_n(dataRowIn, x1, dataRowOut);
+    std::copy_n(dataRowIn, kFeatureSize, dataRowOut + x1);
+    std::copy_n(dataRowIn + x1, input.Width() - x1, dataRowOut + x2);
+    const bool* maskRowIn = mask.ValuePtr(0, y);
+    bool* maskRowOut = mMask.ValuePtr(0, y + kFeatureSize);
+    std::copy_n(maskRowIn, x1, maskRowOut);
+    std::copy_n(maskRowIn, kFeatureSize, maskRowOut + x1);
+    std::copy_n(maskRowIn + x1, input.Width() - x1, maskRowOut + x2);
+    std::fill_n(missing.ValuePtr(x1, y), kFeatureSize, true);
+  }
+}
+
+void removeGap(Mask2D& mask, const Mask2D& mMask) {
+  const size_t x1 = mask.Width() / 2;
+  const size_t x2 = x1 + kFeatureSize;
+  const size_t y1 = mask.Height() / 2;
+  for (size_t y = 0; y != y1; ++y) {
+    const bool* maskRowIn = mMask.ValuePtr(0, y);
+    bool* maskRowOut = mask.ValuePtr(0, y);
+    std::copy_n(maskRowIn, x1, maskRowOut);
+    std::copy_n(maskRowIn + x2, mask.Width() - x1, maskRowOut + x1);
+  }
+  for (size_t y = y1; y != mask.Height(); ++y) {
+    const bool* maskRowIn = mMask.ValuePtr(0, y + kFeatureSize);
+    bool* maskRowOut = mask.ValuePtr(0, y);
+    std::copy_n(maskRowIn, x1, maskRowOut);
+    std::copy_n(maskRowIn + x2, mask.Width() - x1, maskRowOut + x1);
+  }
+}
+
+}  // namespace test_tools
diff --git a/test/algorithms/testtools.h b/test/algorithms/testtools.h
new file mode 100644
index 00000000..0d727abe
--- /dev/null
+++ b/test/algorithms/testtools.h
@@ -0,0 +1,28 @@
+#ifndef TEST_ALGORITHMS_TESTTOOLS_H
+#define TEST_ALGORITHMS_TESTTOOLS_H
+
+#include "../../structures/image2d.h"
+#include "../../structures/mask2d.h"
+
+#include <functional>
+#include <set>
+
+namespace test_tools {
+
+void compareVerticalAlgorithm(
+    std::function<void(const Image2D*, Mask2D*, Mask2D*, size_t, num_t)>
+        executeAlgorithm,
+    const std::set<size_t>& skips);
+
+void compareHorizontalAlgorithm(
+    std::function<void(const Image2D*, Mask2D*, Mask2D*, size_t, num_t)>
+        executeAlgorithm);
+
+void introduceGap(const Image2D& input, const Mask2D& mask, Image2D& mInput,
+                  Mask2D& mMask, Mask2D& missing);
+
+void removeGap(Mask2D& mask, const Mask2D& mMask);
+
+}  // namespace test_tools
+
+#endif
diff --git a/test/experiments/defaultstrategyspeedtest.cpp b/test/experiments/defaultstrategyspeedtest.cpp
index 6d57d610..c7a741ab 100644
--- a/test/experiments/defaultstrategyspeedtest.cpp
+++ b/test/experiments/defaultstrategyspeedtest.cpp
@@ -4,6 +4,7 @@
 #include "../../algorithms/thresholdconfig.h"
 #include "../../algorithms/siroperator.h"
 #include "../../algorithms/sumthreshold.h"
+#include "../../algorithms/sumthresholdmissing.h"
 
 #include "../../lua/telescopefile.h"
 
@@ -17,26 +18,46 @@
 using algorithms::BackgroundTestSet;
 using algorithms::RFITestSet;
 using algorithms::SumThreshold;
+using algorithms::SumThresholdMissing;
 using algorithms::TestSetGenerator;
 using algorithms::ThresholdConfig;
 
+namespace {
+constexpr size_t nRepeats = 3;
+constexpr size_t width = 50000, height = 256;
+}  // namespace
+
 BOOST_AUTO_TEST_SUITE(default_strategy_speed,
                       *boost::unit_test::label("experiment") *
                           boost::unit_test::disabled())
 
 TimeFrequencyData prepareData() {
-  const unsigned width = 50000, height = 256;
-  TimeFrequencyData xx =
-    TestSetGenerator::MakeTestSet(RFITestSet::GaussianBursts, BackgroundTestSet::Empty, width, height);
-  TimeFrequencyData xy =
-    TestSetGenerator::MakeTestSet(RFITestSet::GaussianBursts, BackgroundTestSet::Empty, width, height);
-  TimeFrequencyData yx =
-    TestSetGenerator::MakeTestSet(RFITestSet::GaussianBursts, BackgroundTestSet::Empty, width, height);
-  TimeFrequencyData yy =
-    TestSetGenerator::MakeTestSet(RFITestSet::GaussianBursts, BackgroundTestSet::Empty, width, height);
-  TimeFrequencyData data = TimeFrequencyData::FromLinear(
-      xx.GetImage(0), xx.GetImage(1), xy.GetImage(0), xy.GetImage(1), yx.GetImage(0), yx.GetImage(1), yy.GetImage(0), yy.GetImage(1));
-  return data;
+  TimeFrequencyData data = TestSetGenerator::MakeTestSet(
+      RFITestSet::GaussianBursts, BackgroundTestSet::Empty, width, height);
+  return data.Make(TimeFrequencyData::AmplitudePart);
+}
+
+BOOST_AUTO_TEST_CASE(time_masked_sumthreshold_full) {
+  TimeFrequencyData data = prepareData();
+
+  ThresholdConfig config;
+  config.InitializeLengthsDefault(9);
+  const num_t stddev = data.GetSingleImage()->GetStdDev();
+  config.InitializeThresholdsFromFirstThreshold(6.0 * stddev,
+                                                ThresholdConfig::Rayleigh);
+
+  Image2DCPtr image = data.GetSingleImage();
+  Mask2D zero = Mask2D::MakeSetMask<false>(image->Width(), image->Height());
+  Mask2D mask(zero);
+
+  double t = 0.0;
+  for (size_t i = 0; i != 10; ++i) {
+    Stopwatch watch(true);
+    config.ExecuteWithMissing(image.get(), &mask, &zero, true, 1.0, 1.0);
+    t += watch.Seconds();
+    mask = zero;
+  }
+  std::cout << "10 runs: " << t << '\n';
 }
 
 BOOST_AUTO_TEST_CASE(time_sumthreshold_n_hori) {
@@ -44,11 +65,11 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_hori) {
 
   ThresholdConfig config;
   config.InitializeLengthsDefault(9);
-  num_t stddev = data.GetSingleImage()->GetStdDev();
+  const num_t stddev = data.GetSingleImage()->GetStdDev();
   config.InitializeThresholdsFromFirstThreshold(6.0 * stddev,
                                                 ThresholdConfig::Rayleigh);
-  const size_t N = 100;
-  double hor = 0.0, sseHor = 0.0, avxHorDumas = 0.0, selectedHor = 0.0;
+  double hor = 0.0, sseHor = 0.0, avxHorDumas = 0.0, selectedHor = 0.0,
+         missing = 0.0;
   for (unsigned i = 0; i < 9; ++i) {
     const unsigned length = config.GetHorizontalLength(i);
     const double threshold = config.GetHorizontalThreshold(i);
@@ -57,8 +78,9 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_hori) {
     SumThreshold::VerticalScratch vScratch(input->Width(), input->Height());
 
     Mask2D mask(*data.GetSingleMask()), maskInp;
+    Mask2D zero = Mask2D::MakeSetMask<false>(mask.Width(), mask.Height());
     Stopwatch watch(true);
-    for (size_t j = 0; j != N; ++j) {
+    for (size_t j = 0; j != nRepeats; ++j) {
       maskInp = mask;
       SumThreshold::HorizontalLargeReference(input.get(), &maskInp, &scratch,
                                              length, threshold);
@@ -67,10 +89,21 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_hori) {
     Logger::Info << "Horizontal, length " << length << ": " << watch.ToString()
                  << '\n';
 
+    mask = *data.GetSingleMask();
+    watch.Reset(true);
+    for (size_t j = 0; j != nRepeats; ++j) {
+      maskInp = mask;
+      SumThresholdMissing::Horizontal(*input, maskInp, zero, scratch, length,
+                                      threshold);
+    }
+    missing += watch.Seconds();
+    Logger::Info << "Horizontal with missing, length " << length << ": "
+                 << watch.ToString() << '\n';
+
 #ifdef __SSE__
     mask = *data.GetSingleMask();
     watch.Reset(true);
-    for (size_t j = 0; j != N; ++j) {
+    for (size_t j = 0; j != nRepeats; ++j) {
       maskInp = mask;
       SumThreshold::HorizontalLargeSSE(input.get(), &maskInp, &scratch, length,
                                        threshold);
@@ -83,7 +116,7 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_hori) {
 #ifdef __AVX2__
     mask = *data.GetSingleMask();
     watch.Reset(true);
-    for (size_t j = 0; j != N; ++j) {
+    for (size_t j = 0; j != nRepeats; ++j) {
       maskInp = mask;
       SumThreshold::HorizontalAVXDumas(input.get(), &maskInp, length,
                                        threshold);
@@ -95,7 +128,7 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_hori) {
 
     mask = *data.GetSingleMask();
     watch.Reset(true);
-    for (size_t j = 0; j != N; ++j) {
+    for (size_t j = 0; j != nRepeats; ++j) {
       maskInp = mask;
       SumThreshold::HorizontalLarge(input.get(), &maskInp, &scratch, length,
                                     threshold);
@@ -106,6 +139,7 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_hori) {
 
     Logger::Info << "Summed values:\n"
                  << "- Horizontal ref  : " << hor << "\n"
+                 << "- Horizontal missing  : " << missing << "\n"
                  << "- Horizontal SSE  : " << sseHor << "\n"
                  << "- Horizontal AVX D: " << avxHorDumas << "\n"
                  << "- Selected horiz  : " << selectedHor << "\n";
@@ -117,11 +151,11 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_vert) {
 
   ThresholdConfig config;
   config.InitializeLengthsDefault(9);
-  num_t stddev = data.GetSingleImage()->GetStdDev();
+  const num_t stddev = data.GetSingleImage()->GetStdDev();
   config.InitializeThresholdsFromFirstThreshold(6.0 * stddev,
                                                 ThresholdConfig::Rayleigh);
-  const size_t N = 100;
   double vert = 0.0, sseVert = 0.0, avxVert = 0.0, avxVertDumas = 0.0,
+         missingRef = 0.0, missingStacked = 0.0, missingConsecutive = 0.0,
          selectedVer = 0.0;
   for (unsigned i = 0; i < 9; ++i) {
     const unsigned length = config.GetHorizontalLength(i);
@@ -131,9 +165,10 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_vert) {
     SumThreshold::VerticalScratch vScratch(input->Width(), input->Height());
 
     Mask2D mask(*data.GetSingleMask()), maskInp;
+    Mask2D zero = Mask2D::MakeSetMask<false>(mask.Width(), mask.Height());
 
     Stopwatch watch(true);
-    for (size_t j = 0; j != N; ++j) {
+    for (size_t j = 0; j != nRepeats; ++j) {
       maskInp = mask;
       SumThreshold::VerticalLargeReference(input.get(), &maskInp, &scratch,
                                            length, threshold);
@@ -142,10 +177,45 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_vert) {
     Logger::Info << "Vertical, length " << length << ": " << watch.ToString()
                  << '\n';
 
+    mask = *data.GetSingleMask();
+    watch.Reset(true);
+    for (size_t j = 0; j != nRepeats; ++j) {
+      maskInp = mask;
+      SumThresholdMissing::VerticalReference(*input, maskInp, zero, scratch,
+                                             length, threshold);
+    }
+    missingRef += watch.Seconds();
+    Logger::Info << "Vertical missing ref, length " << length << ": "
+                 << watch.ToString() << '\n';
+
+    mask = *data.GetSingleMask();
+    watch.Reset(true);
+    SumThresholdMissing::VerticalCache vCache;
+    SumThresholdMissing::InitializeVertical(vCache, *input, zero);
+    for (size_t j = 0; j != nRepeats; ++j) {
+      maskInp = mask;
+      SumThresholdMissing::VerticalStacked(vCache, *input, maskInp, zero,
+                                           scratch, length, threshold);
+    }
+    missingStacked += watch.Seconds();
+    Logger::Info << "Vertical missing stacked, length " << length << ": "
+                 << watch.ToString() << '\n';
+
+    mask = *data.GetSingleMask();
+    watch.Reset(true);
+    for (size_t j = 0; j != nRepeats; ++j) {
+      maskInp = mask;
+      SumThresholdMissing::VerticalConsecutive(*input, maskInp, zero, scratch,
+                                               length, threshold);
+    }
+    missingConsecutive += watch.Seconds();
+    Logger::Info << "Vertical missing consecutive, length " << length << ": "
+                 << watch.ToString() << '\n';
+
 #ifdef __SSE__
     mask = *data.GetSingleMask();
     watch.Reset(true);
-    for (size_t j = 0; j != N; ++j) {
+    for (size_t j = 0; j != nRepeats; ++j) {
       maskInp = mask;
       SumThreshold::VerticalLargeSSE(input.get(), &maskInp, &scratch, length,
                                      threshold);
@@ -158,7 +228,7 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_vert) {
 #ifdef __AVX2__
     mask = *data.GetSingleMask();
     watch.Reset(true);
-    for (size_t j = 0; j != N; ++j) {
+    for (size_t j = 0; j != nRepeats; ++j) {
       maskInp = mask;
       SumThreshold::VerticalLargeAVX(input.get(), &maskInp, &scratch, length,
                                      threshold);
@@ -169,7 +239,7 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_vert) {
 
     mask = *data.GetSingleMask();
     watch.Reset(true);
-    for (size_t j = 0; j != N; ++j) {
+    for (size_t j = 0; j != nRepeats; ++j) {
       maskInp = mask;
       SumThreshold::VerticalAVXDumas(input.get(), &maskInp, &vScratch, length,
                                      threshold);
@@ -181,7 +251,7 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_vert) {
 
     mask = *data.GetSingleMask();
     watch.Reset(true);
-    for (size_t j = 0; j != N; ++j) {
+    for (size_t j = 0; j != nRepeats; ++j) {
       maskInp = mask;
       SumThreshold::VerticalLarge(input.get(), &maskInp, &scratch, &vScratch,
                                   length, threshold);
@@ -191,11 +261,14 @@ BOOST_AUTO_TEST_CASE(time_sumthreshold_n_vert) {
                  << watch.ToString() << '\n';
 
     Logger::Info << "Summed values:\n"
-                 << "- Vertical ref    : " << vert << "\n"
-                 << "- Vertical SSE    : " << sseVert << "\n"
-                 << "- Vertical AVX    : " << avxVert << "\n"
-                 << "- Vertical AVX D  : " << avxVertDumas << "\n"
-                 << "- Selected vertic : " << selectedVer << "\n";
+                 << "- Vertical ref        : " << vert << "\n"
+                 << "- Vertical missing ref: " << missingRef << "\n"
+                 << "- Vertical m. stacked : " << missingStacked << "\n"
+                 << "- Vertical m. consec  : " << missingConsecutive << "\n"
+                 << "- Vertical SSE        : " << sseVert << "\n"
+                 << "- Vertical AVX        : " << avxVert << "\n"
+                 << "- Vertical AVX D      : " << avxVertDumas << "\n"
+                 << "- Selected vertic     : " << selectedVer << "\n";
   }
 }
 
diff --git a/test/interface/interfacetest.cpp b/test/interface/interfacetest.cpp
index e9cf62e0..909f70e7 100644
--- a/test/interface/interfacetest.cpp
+++ b/test/interface/interfacetest.cpp
@@ -59,7 +59,9 @@ BOOST_AUTO_TEST_CASE(run_default_strategy_with_input) {
 BOOST_AUTO_TEST_CASE(runs) {
   size_t width = 200, height = 50;
   Mask2D gtMask = Mask2D::MakeSetMask<false>(width, height);
-  TimeFrequencyData data = TestSetGenerator::MakeTestSet(algorithms::RFITestSet::FullBandBursts, algorithms::BackgroundTestSet::Empty, width, height);
+  TimeFrequencyData data = TestSetGenerator::MakeTestSet(
+      algorithms::RFITestSet::FullBandBursts,
+      algorithms::BackgroundTestSet::Empty, width, height);
   Image2DCPtr real = data.GetImage(0);
   Image2DCPtr imag = data.GetImage(1);
   aoflagger::AOFlagger flagger;
diff --git a/test/lua/defaultstrategytest.cpp b/test/lua/defaultstrategytest.cpp
index 3d3523c8..d44f581c 100644
--- a/test/lua/defaultstrategytest.cpp
+++ b/test/lua/defaultstrategytest.cpp
@@ -51,19 +51,20 @@ static std::pair<double, double> checkResult(TimeFrequencyData& data,
 
 BOOST_AUTO_TEST_CASE(stokesi_complex) {
   size_t width = 200, height = 50;
-  TimeFrequencyData data =
-    TestSetGenerator::MakeTestSet(RFITestSet::FullBandBursts, BackgroundTestSet::Empty, width, height);
+  TimeFrequencyData data = TestSetGenerator::MakeTestSet(
+      RFITestSet::FullBandBursts, BackgroundTestSet::Empty, width, height);
 
   auto fpfn = checkResult(data, *data.GetSingleMask());
-  BOOST_CHECK_LT(fpfn.first, 5e-3);   // False positives-ratio should be < 0.5 %
-  BOOST_CHECK_LT(fpfn.second, 0.015);  // False negatives-ratio should be < 1.5 %
+  BOOST_CHECK_LT(fpfn.first, 5e-3);  // False positives-ratio should be < 0.5 %
+  BOOST_CHECK_LT(fpfn.second,
+                 0.015);  // False negatives-ratio should be < 1.5 %
 }
 
 BOOST_AUTO_TEST_CASE(stokesi_amplitude) {
   size_t width = 200, height = 50;
   Mask2D mask = Mask2D::MakeSetMask<false>(width, height);
-  TimeFrequencyData data =
-    TestSetGenerator::MakeTestSet(RFITestSet::FullBandBursts, BackgroundTestSet::Empty, width, height);
+  TimeFrequencyData data = TestSetGenerator::MakeTestSet(
+      RFITestSet::FullBandBursts, BackgroundTestSet::Empty, width, height);
   Image2DCPtr amplitude = data.GetImage(0);
 
   TimeFrequencyData tfData(TimeFrequencyData::AmplitudePart,
@@ -80,12 +81,11 @@ BOOST_AUTO_TEST_CASE(full_polarization) {
   PolarizationEnum pols[4] = {Polarization::XX, Polarization::XY,
                               Polarization::YX, Polarization::YY};
   for (size_t i = 0; i != 4; ++i) {
-    TimeFrequencyData real =
-      TestSetGenerator::MakeTestSet(RFITestSet::FullBandBursts, BackgroundTestSet::Empty, width, height);
-    TimeFrequencyData imag =
-      TestSetGenerator::MakeTestSet(RFITestSet::FullBandBursts, BackgroundTestSet::Empty, width, height);
-    TimeFrequencyData tfData(pols[i], real.GetImage(0),
-                             imag.GetImage(0));
+    TimeFrequencyData real = TestSetGenerator::MakeTestSet(
+        RFITestSet::FullBandBursts, BackgroundTestSet::Empty, width, height);
+    TimeFrequencyData imag = TestSetGenerator::MakeTestSet(
+        RFITestSet::FullBandBursts, BackgroundTestSet::Empty, width, height);
+    TimeFrequencyData tfData(pols[i], real.GetImage(0), imag.GetImage(0));
     tfData.SetGlobalMask(real.GetMask(0));
     if (i == 0)
       allData = tfData;
-- 
GitLab