brute_force_pairwise.h

Go to the documentation of this file.

#pragma once
 
#include <cstddef>
#include <cstdint>
#include <type_traits>
#include <vector>
 
#include "clustering/math/pairwise.h"
#include "clustering/math/thread.h"
#include "clustering/ndarray.h"
 
namespace clustering {

template <class T> class BruteForcePairwise {
public:
  static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>,
                "BruteForcePairwise<T> requires T to be float or double");

  explicit BruteForcePairwise(const NDArray<T, 2> &points) noexcept : m_points(points) {}

  [[nodiscard]] std::vector<std::vector<std::int32_t>> query(T radius, math::Pool pool) const {
    const std::size_t n = m_points.dim(0);
    std::vector<std::vector<std::int32_t>> adj(n);
    if (n == 0) {
      return adj;
    }
 
    // Reserve a small floor per row so the first push_backs do not trigger the vector-doubling
    // reallocation cascade that otherwise dominates adjacency construction on dense fixtures.
    for (auto &v : adj) {
      v.reserve(16);
    }
 
    const T radiusSq = radius * radius;
    // Symmetric eps-neighbour graph: the kernel emits each unique upper-triangular cell once
    // with `row <= col`. The kernel-side emit only touches `adj[row]` so workers writing
    // disjoint row chunks remain race-free; the mirror push (`adj[col]`.push(row)) runs as a
    // single-threaded post-pass below. Halves the pairwise compute on the brute-force path.
    auto emit = [&adj](std::size_t row, std::size_t col) {
      adj[row].push_back(static_cast<std::int32_t>(col));
    };
    math::pairwiseSqEuclideanThresholdedSymmetric(m_points, radiusSq, pool, emit);
 
    // Mirror pass: each surviving upper-triangular pair `(i, j)` lives in `adj[i]`; push
    // @c i to `adj[j]` sequentially so the parallel sweep above never crosses worker
    // boundaries. The size snapshot stops the loop from walking the freshly-mirrored entries
    // when the next outer @c i lands on a row that earlier iterations already augmented.
    for (std::size_t i = 0; i < n; ++i) {
      const std::size_t origSize = adj[i].size();
      for (std::size_t k = 0; k < origSize; ++k) {
        const auto jSigned = adj[i][k];
        const auto j = static_cast<std::size_t>(jSigned);
        if (j > i) {
          adj[j].push_back(static_cast<std::int32_t>(i));
        }
      }
    }
    return adj;
  }
 
private:
  const NDArray<T, 2> &m_points; 
};
 
} // namespace clustering