brute_force_pairwise.h

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <type_traits>
#include <utility>
#include <vector>
 
#include "clustering/index/range_query.h"
#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]] index::CoreAdjacency query(T radius, std::size_t minPts, math::Pool pool) const {
    const std::size_t n = m_points.dim(0);
    index::CoreAdjacency out;
    out.rows.resize(n);
    out.isCore.assign(n, 0);
    if (n == 0) {
      return out;
    }
    std::vector<std::vector<std::int32_t>> &adj = out.rows;
 
    const std::size_t d = m_points.dim(1);
    std::size_t adjReserveFloor = 16;
    if (d == 32) {
      adjReserveFloor = 24;
    } else if (d == 64) {
      adjReserveFloor = 20;
    }
    // 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.
    // Fanning the reserves out matters at high worker counts: the row allocations are a
    // serial malloc train that idles every worker before the sweep starts, and the allocator's
    // per-thread arenas let the fan-out scale.
    const auto reserveRows = [&](std::size_t lo, std::size_t hi) {
      for (std::size_t i = lo; i < hi; ++i) {
        adj[i].reserve(adjReserveFloor);
      }
    };
    const bool fanOut = pool.shouldParallelize(n, 256, 2);
    if (fanOut) {
      pool.parallelForBlocks(std::size_t{0}, n, std::size_t{0}, reserveRows);
    } else {
      reserveRows(0, n);
    }
 
    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.
    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 degrees without mirror edges: workers accumulate `j < i` neighbour counts in
    // private histograms that stay cache-resident, then a row-partitioned reduction folds them.
    // Scattering the mirrored edges themselves would be a latency-bound write per edge into a
    // random destination row; the degree is all the core verdict needs.
    const std::size_t workers = pool.workerCount();
    std::vector<std::vector<std::uint32_t>> workerCounts(workers);
    const auto countRange = [&](std::size_t lo, std::size_t hi) {
      std::vector<std::uint32_t> &counts = workerCounts[math::Pool::workerIndex()];
      if (counts.empty()) {
        counts.assign(n, 0);
      }
      for (std::size_t i = lo; i < hi; ++i) {
        for (const std::int32_t neighbor : adj[i]) {
          const auto j = static_cast<std::size_t>(neighbor);
          if (j > i) {
            ++counts[j];
          }
        }
      }
    };
    std::vector<std::uint32_t> mirrorDeg(n, 0);
    const auto reduceRange = [&](std::size_t lo, std::size_t hi) {
      for (const std::vector<std::uint32_t> &counts : workerCounts) {
        if (counts.empty()) {
          continue;
        }
        for (std::size_t i = lo; i < hi; ++i) {
          mirrorDeg[i] += counts[i];
        }
      }
    };
    const auto flagRange = [&](std::size_t lo, std::size_t hi) {
      for (std::size_t i = lo; i < hi; ++i) {
        out.isCore[i] =
            (adj[i].size() + mirrorDeg[i] >= minPts) ? std::uint8_t{1} : std::uint8_t{0};
      }
    };
    if (fanOut) {
      pool.parallelForBlocks(std::size_t{0}, n, std::size_t{0}, countRange);
      pool.parallelForBlocks(std::size_t{0}, n, std::size_t{0}, reduceRange);
      pool.parallelForBlocks(std::size_t{0}, n, std::size_t{0}, flagRange);
    } else {
      countRange(0, n);
      reduceRange(0, n);
      flagRange(0, n);
    }
 
    // Complete the non-core rows only. Workers collect the surviving `(dest, src)` mirrors
    // into private lists; the survivor set is a vanishing fraction of the edge set, so the
    // ordered serial apply costs a handful of pushes and keeps row contents deterministic.
    std::vector<std::vector<std::pair<std::int32_t, std::int32_t>>> workerMirrors(workers);
    const auto collectRange = [&](std::size_t lo, std::size_t hi) {
      std::vector<std::pair<std::int32_t, std::int32_t>> &mirrors =
          workerMirrors[math::Pool::workerIndex()];
      for (std::size_t i = lo; i < hi; ++i) {
        for (const std::int32_t neighbor : adj[i]) {
          const auto j = static_cast<std::size_t>(neighbor);
          if (j > i && out.isCore[j] == 0) {
            mirrors.emplace_back(static_cast<std::int32_t>(j), static_cast<std::int32_t>(i));
          }
        }
      }
    };
    if (fanOut) {
      pool.parallelForBlocks(std::size_t{0}, n, std::size_t{0}, collectRange);
    } else {
      collectRange(0, n);
    }
    std::vector<std::pair<std::int32_t, std::int32_t>> mirrors;
    for (const auto &local : workerMirrors) {
      mirrors.insert(mirrors.end(), local.begin(), local.end());
    }
    std::sort(mirrors.begin(), mirrors.end());
    for (const auto &[dest, src] : mirrors) {
      adj[static_cast<std::size_t>(dest)].push_back(src);
    }
    return out;
  }
 
private:
  const NDArray<T, 2> &m_points; 
};
 
} // namespace clustering