Handle histogram edge cases on ARM hardware

Author: jngrad

.github/workflows/push_pull.yml

@@ -40,6 +40,9 @@ jobs:
           build_procs: 3
           check_procs: 3
           with_ccache: 'true'
+          with_walberla: 'true'
+          with_walberla_fft: 'false'
+          with_walberla_avx: 'false'
 
   debian:
     runs-on: ubuntu-latest

maintainer/CI/build_cmake.sh

@@ -132,6 +132,7 @@ set_default_value with_fftw true
 set_default_value with_gsl true
 set_default_value with_scafacos false
 set_default_value with_walberla false
+set_default_value with_walberla_fft true
 set_default_value with_walberla_avx false
 set_default_value with_stokesian_dynamics false
 set_default_value test_timeout 500
@@ -168,7 +169,9 @@ cmake_params="${cmake_params} -D ESPRESSO_BUILD_WITH_STOKESIAN_DYNAMICS=${with_s
 cmake_params="${cmake_params} -D ESPRESSO_BUILD_WITH_WALBERLA=${with_walberla}"
 
 if [ "${with_walberla}" = true ]; then
-  cmake_params="${cmake_params} -D ESPRESSO_BUILD_WITH_WALBERLA_FFT=ON"
+  if [ "${with_walberla_fft}" = true ]; then
+    cmake_params="${cmake_params} -D ESPRESSO_BUILD_WITH_WALBERLA_FFT=ON"
+  fi
   if [ "${with_walberla_avx}" = true ]; then
     cmake_params="${cmake_params} -D ESPRESSO_BUILD_WITH_WALBERLA_AVX=ON"
   fi

src/utils/include/utils/Histogram.hpp

@@ -105,13 +105,10 @@ class Histogram {
       throw std::invalid_argument("Wrong dimensions for the value");
     }
     if (check_limits(pos)) {
-      boost::array<array_index, M + 1> index;
-      for (std::size_t i = 0; i < M; ++i) {
-        index[i] = calc_bin_index(pos[i], m_limits[i].first, m_bin_sizes[i]);
-      }
-      for (array_index i = 0; i < static_cast<array_index>(N); ++i) {
-        index.back() = i;
-        m_array(index) += value[static_cast<std::size_t>(i)];
+      auto index = calc_bin_index(pos);
+      for (std::size_t i = 0; i < N; ++i) {
+        index.back() = static_cast<array_index>(i);
+        m_array(index) += value[i];
         m_count(index)++;
       }
     }
@@ -121,20 +118,42 @@ class Histogram {
   virtual void normalize() {
     auto const bin_volume = std::accumulate(
         m_bin_sizes.begin(), m_bin_sizes.end(), U{1}, std::multiplies<U>());
-    std::transform(
-        m_array.data(), m_array.data() + m_array.num_elements(), m_array.data(),
+    std::ranges::transform(
+        std::span(m_array.data(), m_array.num_elements()), m_array.data(),
         [bin_volume](T v) { return static_cast<T>(v / bin_volume); });
   }
 
+protected:
+  // internal function exposed for unit testing purposes
+  long correct_off_by_one(long bin, long n_bins) const {
+    if (bin == -1l) {
+      ++bin;
+    } else if (bin == n_bins) {
+      --bin;
+    }
+    return bin;
+  }
+
 private:
   /**
    * \brief Calculate the bin index.
-   * \param value  Position on that dimension.
-   * \param offset Bin offset on that dimension.
-   * \param size   Bin size on that dimension.
+   * \param pos  Position.
    */
-  array_index calc_bin_index(double value, double offset, double size) const {
-    return static_cast<array_index>(std::floor((value - offset) / size));
+  auto calc_bin_index(std::span<const U> const &pos) const {
+    boost::array<array_index, M + 1> index;
+    for (std::size_t i = 0; i < M; ++i) {
+      auto const offset = m_limits[i].first;
+      auto const size = m_bin_sizes[i];
+      auto const n_bins = static_cast<long>(m_n_bins[i]);
+      auto const bin = static_cast<long>(std::floor((pos[i] - offset) / size));
+      // handle edge cases when the position is exactly between two bins:
+      // due to precision loss in the offset subtraction, the bin index might
+      // be off by one, so we fold it here back inside the valid range
+      auto const folded_bin = correct_off_by_one(bin, n_bins);
+      assert((folded_bin == bin) or check_limits(pos));
+      index[i] = static_cast<array_index>(folded_bin);
+    }
+    return index;
   }
 
   /**
@@ -153,7 +172,7 @@ class Histogram {
    * \brief Check if the position lies within the histogram limits.
    * \param pos     Position to check.
    */
-  bool check_limits(std::span<const U> pos) const {
+  bool check_limits(std::span<const U> const &pos) const {
     assert(pos.size() == M);
     bool within_range = true;
     for (std::size_t i = 0; i < M; ++i) {
@@ -165,7 +184,7 @@ class Histogram {
 
   std::array<std::size_t, M + 1> m_array_dim() const {
     std::array<std::size_t, M + 1> dimensions;
-    std::copy(m_n_bins.begin(), m_n_bins.end(), dimensions.begin());
+    std::ranges::copy(m_n_bins, dimensions.begin());
     dimensions.back() = N;
     return dimensions;
   }
@@ -193,14 +212,15 @@ class Histogram {
  */
 template <typename T, std::size_t N, std::size_t M = 3, typename U = double>
 class CylindricalHistogram : public Histogram<T, N, M, U> {
-  using Histogram<T, N, M, U>::m_n_bins;
-  using Histogram<T, N, M, U>::m_limits;
-  using Histogram<T, N, M, U>::m_bin_sizes;
-  using Histogram<T, N, M, U>::m_array;
-  using typename Histogram<T, N, M, U>::array_index;
+  using Base = Histogram<T, N, M, U>;
+  using Base::m_array;
+  using Base::m_bin_sizes;
+  using Base::m_limits;
+  using Base::m_n_bins;
+  using typename Base::array_index;
 
 public:
-  using Histogram<T, N, M, U>::Histogram;
+  using Base::Histogram;
 
   void normalize() override {
     auto const min_r = m_limits[0].first;
@@ -214,8 +234,8 @@ class CylindricalHistogram : public Histogram<T, N, M, U> {
       auto const bin_volume = (r_right * r_right - r_left * r_left) *
                               z_bin_size * phi_bin_size / U(2);
       auto *begin = m_array[i].origin();
-      std::transform(
-          begin, begin + m_array[i].num_elements(), begin,
+      std::ranges::transform(
+          std::span(begin, m_array[i].num_elements()), begin,
           [bin_volume](T v) { return static_cast<T>(v / bin_volume); });
     }
   }

src/utils/tests/histogram_test.cpp

@@ -105,3 +105,45 @@ BOOST_AUTO_TEST_CASE(cylindrical_histogram) {
                                 std::vector<double>{{0.0, 0.0}}),
                     std::invalid_argument);
 }
+
+BOOST_AUTO_TEST_CASE(histogram_bin_folding) {
+  class TestHistogram : public Utils::Histogram<double, 2, 2, double> {
+    using Base = Histogram<double, 2, 2, double>;
+
+  public:
+    using Base::correct_off_by_one;
+    using Base::Histogram;
+  };
+  std::array<std::size_t, 2> n_bins{{10, 10}};
+  std::array<std::pair<double, double>, 2> limits{
+      {std::make_pair(1.0, 20.0), std::make_pair(5.0, 10.0)}};
+  auto hist = TestHistogram(n_bins, limits);
+  auto old_bin = 0l;
+  auto new_bin = 0l;
+  auto constexpr max_index = 2l;
+  auto constexpr min_index = 0l;
+  // no folding from the right
+  old_bin = max_index - 1l;
+  new_bin = hist.correct_off_by_one(old_bin, max_index);
+  BOOST_CHECK_EQUAL(new_bin, old_bin);
+  // folding from the right
+  old_bin = max_index;
+  new_bin = hist.correct_off_by_one(old_bin, max_index);
+  BOOST_CHECK_EQUAL(new_bin, old_bin - 1l);
+  // no folding from the right
+  old_bin = max_index + 1l;
+  new_bin = hist.correct_off_by_one(old_bin, max_index);
+  BOOST_CHECK_EQUAL(new_bin, old_bin);
+  // no folding from the left
+  old_bin = min_index;
+  new_bin = hist.correct_off_by_one(old_bin, max_index);
+  BOOST_CHECK_EQUAL(new_bin, old_bin);
+  // folding from the left
+  old_bin = min_index - 1l;
+  new_bin = hist.correct_off_by_one(old_bin, max_index);
+  BOOST_CHECK_EQUAL(new_bin, old_bin + 1l);
+  // no folding from the left
+  old_bin = min_index - 2l;
+  new_bin = hist.correct_off_by_one(old_bin, max_index);
+  BOOST_CHECK_EQUAL(new_bin, old_bin);
+}