SortingThresholdBenchmark.h

Go to the documentation of this file.

 
#pragma once
 
#include <array>
#include <cmath>
#include <memory>
#include <random>
 
#include "autopas/baseFunctors/CellFunctor.h"
#include "autopas/options/SortingDirectionOption.h"
#include "autopas/utils/SortingThresholdInfo2B.h"
#include "autopas/utils/Timer.h"
#include "autopas/utils/generators/UniformGenerator.h"
#include "autopas/utils/logging/Logger.h"
 
namespace autopas {
 
class SortingThresholdBenchmark {
 public:
  SortingThresholdBenchmark(size_t aosSortingThresholdDefault, size_t soaSortingThresholdDefault)
      : _soaThresholds(std::make_shared<const SortingThresholdInfoSingle>(soaSortingThresholdDefault)),
        _aosThresholds(std::make_shared<const SortingThresholdInfoSingle>(aosSortingThresholdDefault)) {}
 
  [[nodiscard]] std::shared_ptr<const SortingThresholdInfoInterface> getSoAThreshold() const { return _soaThresholds; }
 
  [[nodiscard]] std::shared_ptr<const SortingThresholdInfoInterface> getAoSThreshold() const { return _aosThresholds; }
 
  [[nodiscard]] bool hasRunSoA() const { return _hasRunSoA; }
 
  [[nodiscard]] bool hasRunAoS() const { return _hasRunAoS; }
 
  template <class Functor_T, class Particle_T>
  void runSoABenchmark(Functor_T &functor, const Particle_T &defaultParticle) {
    SortingThresholdInfo2B thresholds = SortingThresholdInfo2B(0);
    for (const auto &n3 : Newton3Option::getAllOptions()) {
      for (const auto &cellDirection : SortingDirectionOption::getAllOptions()) {
        thresholds.setThresholdByOption(
            n3, cellDirection, /*bidirectional=*/true,
            runSearch<Functor_T, Particle_T, true, true>(functor, defaultParticle, cellDirection, n3));
        if (n3 == Newton3Option::disabled) {
          thresholds.setThresholdByOption(
              n3, cellDirection, /*bidirectional=*/false,
              runSearch<Functor_T, Particle_T, true, false>(functor, defaultParticle, cellDirection, n3));
        }
      }
    }
    _soaThresholds = std::make_shared<const SortingThresholdInfo2B>(thresholds);
    _hasRunSoA = true;
  }
 
  template <class Functor_T, class Particle_T>
  void runAoSBenchmark(Functor_T &functor, const Particle_T &defaultParticle) {
    SortingThresholdInfo2B thresholds = SortingThresholdInfo2B(0);
    for (const auto &n3 : Newton3Option::getAllOptions()) {
      for (const auto &cellDirection : SortingDirectionOption::getAllOptions()) {
        thresholds.setThresholdByOption(
            n3, cellDirection, /*bidirectional=*/true,
            runSearch<Functor_T, Particle_T, false, true>(functor, defaultParticle, cellDirection, n3));
        if (n3 == Newton3Option::disabled) {
          thresholds.setThresholdByOption(
              n3, cellDirection, /*bidirectional=*/false,
              runSearch<Functor_T, Particle_T, false, false>(functor, defaultParticle, cellDirection, n3));
        }
      }
    }
    _aosThresholds = std::make_shared<const SortingThresholdInfo2B>(thresholds);
    _hasRunAoS = true;
  }
 
 private:
  bool _hasRunSoA{false};
 
  bool _hasRunAoS{false};
 
  std::shared_ptr<const SortingThresholdInfoInterface> _soaThresholds;
 
  std::shared_ptr<const SortingThresholdInfoInterface> _aosThresholds;
 
  const size_t _iterations = 25;
 
  const size_t _repetitions = 100;
  const size_t _maxSoAParticles = 250;
 
  const size_t _maxAoSParticles = 50;
 
  const double _sortedWinMarginFraction = 0.05;
 
  const double _requiredSortedWinRatio = 0.7;
 
  const double _scatterFactor = 0.2;
 
  template <class Functor_T, class Particle_T, bool useSoA, bool bidirectional>
  size_t executeRun(Functor_T &functor, const Particle_T &defaultParticle, SortingDirectionOption cellDirection,
                    size_t numParticles, Newton3Option newton3) {
    using BenchCell = FullParticleCell<Particle_T>;
    using BenchCF = internal::CellFunctor<BenchCell, Functor_T, bidirectional>;
 
    const double cutoff = functor.getCutoff();
    const double invSqrt3 = 1. / sqrt(3.);
    const double invSqrt2 = 1. / sqrt(2.);
 
    auto numParticlesEqDistr = static_cast<size_t>(std::ceil(numParticles * (1. - _scatterFactor)));
    size_t numParticlesScatter = numParticles - numParticlesEqDistr + numParticlesEqDistr % 2;
    size_t numParticlesPerCell = numParticlesEqDistr / 2;
 
    // Prepare for random scattering.
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_int_distribution<size_t> distrib(0, numParticlesScatter);
 
    // Setup Cell Layout
    BenchCell cell1, cell2;
 
    std::array cell1Low = {0., 0., 0.};
    std::array cell1High = {cutoff, cutoff, cutoff};
 
    std::array cell2Low = {0., 0., 0.};
    std::array cell2High = {cutoff, cutoff, cutoff};
 
    std::array sortingDirection = {0., 0., 0.};
 
    switch (cellDirection) {
      case SortingDirectionOption::corner:
        cell2Low = {cutoff, cutoff, cutoff};
        cell2High = {2. * cutoff, 2. * cutoff, 2. * cutoff};
        sortingDirection = {invSqrt3, invSqrt3, invSqrt3};
        break;
      case SortingDirectionOption::edge:
        cell2Low = {cutoff, cutoff, 0.};
        cell2High = {2. * cutoff, 2. * cutoff, cutoff};
        sortingDirection = {invSqrt2, invSqrt2, 0.};
        break;
      case SortingDirectionOption::face:
        cell2Low = {cutoff, 0., 0.};
        cell2High = {2 * cutoff, cutoff, cutoff};
        sortingDirection = {1, 0., 0.};
        break;
      default:
        utils::ExceptionHandler::exception("Cell Direction {} is not a valid/supported Direction!", cellDirection);
    }
 
    utils::Timer sortedTimer, unsortedTimer;
    BenchCF cellFunctor{functor, cutoff, useSoA ? DataLayoutOption::soa : DataLayoutOption::aos,
                        newton3 == Newton3Option::enabled};
    size_t sortedWins = 0;
    // Set to 0 so whether sorting happens is controlled entirely through the sorting direction.
    const SortingThresholdInfoSingle zeroThreshold(0);
    cellFunctor.setSoASortingThresholds(zeroThreshold);
    cellFunctor.setAoSSortingThresholds(zeroThreshold);
 
    // For SoA Forces won't accumulate, as SoAExtractor is never called, for AoS they will as AoS modifies the values
    // directly.
    auto measureUnsorted = [&]() {
      unsortedTimer.start();
      for (size_t j = 0; j < _iterations; j++) {
        if constexpr (useSoA) {
          functor.SoALoader(cell1, cell1._particleSoABuffer, 0, false);
          functor.SoALoader(cell2, cell2._particleSoABuffer, 0, false);
        }
        // A sorting direction of (0, 0, 0) disables sorting.
        cellFunctor.processCellPair(cell1, cell2, {0., 0., 0.});
      }
      return unsortedTimer.stop();
    };
    auto measureSorted = [&]() {
      sortedTimer.start();
      for (size_t j = 0; j < _iterations; j++) {
        if constexpr (useSoA) {
          functor.SoALoader(cell1, cell1._particleSoABuffer, 0, false);
          functor.SoALoader(cell2, cell2._particleSoABuffer, 0, false);
        }
        cellFunctor.processCellPair(cell1, cell2, sortingDirection);
      }
      return sortedTimer.stop();
    };
 
    for (size_t i = 0; i < _repetitions; i++) {
      cell1.clear();
      cell2.clear();
 
      // Generate the number of scattered particles for cell1 and cell2
      size_t toAddCell1 = distrib(gen);
      // Vary the seed per repetition per cell so each repetition samples a fresh particle
      // layout instead of repeatedly timing the exact same configuration.
      generators::UniformGenerator::fillWithParticles(cell1, defaultParticle, cell1Low, cell1High,
                                                      numParticlesPerCell + toAddCell1,
                                                      static_cast<unsigned int>(2 * i));
      generators::UniformGenerator::fillWithParticles(cell2, defaultParticle, cell2Low, cell2High,
                                                      numParticlesPerCell + numParticlesScatter - toAddCell1,
                                                      static_cast<unsigned int>(2 * i + 1));
 
      long unsortedDelta = 0;
      long sortedDelta = 0;
      // Alternate which path is measured first. Whichever path runs second inherits warm caches and a settled
      // branch predictor from the first, which would otherwise make it look systematically faster than it is.
      if (i % 2 == 0) {
        unsortedDelta = measureUnsorted();
        sortedDelta = measureSorted();
      } else {
        sortedDelta = measureSorted();
        unsortedDelta = measureUnsorted();
      }
 
      AutoPasLog(TRACE, "SortingThresholdBenchmark rep {}/{} cell direction={} n={}: unsorted={}ns sorted={}ns", i + 1,
                 _repetitions, cellDirection, numParticles, unsortedDelta, sortedDelta);
 
      // A repetition only counts as a "sorted win" if it clears the margin: see _sortedWinMarginFraction.
      if (static_cast<double>(sortedDelta) < static_cast<double>(unsortedDelta) * (1. - _sortedWinMarginFraction)) {
        ++sortedWins;
      }
    }
 
    const long meanSorted = sortedTimer.getTotalTime() / static_cast<long>(_repetitions);
    const long meanUnsorted = unsortedTimer.getTotalTime() / static_cast<long>(_repetitions);
    AutoPasLog(TRACE,
               "SortingThresholdBenchmark cell direction={} n={}: mean unsorted={}ns mean sorted={}ns sortedWins={}/{}",
               cellDirection, numParticles, meanUnsorted, meanSorted, sortedWins, _repetitions);
    return sortedWins;
  }
 
  template <class Functor_T, class Particle_T, bool useSoA, bool bidirectional>
  size_t runSearch(Functor_T &functor, const Particle_T &defaultParticle, SortingDirectionOption cellDirection,
                   Newton3Option newton3) {
    size_t lowCount = 0;
    const size_t upperBound = useSoA ? _maxSoAParticles : _maxAoSParticles;
    size_t highCount = upperBound;
 
    while (lowCount < highCount) {
      size_t mid = lowCount + (highCount - lowCount) / 2;
 
      const auto outcome = executeRun<Functor_T, Particle_T, useSoA, bidirectional>(functor, defaultParticle,
                                                                                    cellDirection, mid, newton3);
      const double winRatio = static_cast<double>(outcome) / static_cast<double>(_repetitions);
 
      // Conservative decision rule: only accept "sorted wins" once a clear majority of repetitions
      // agree by a clear margin (see _sortedWinMarginFraction and _requiredSortedWinRatio).
      if (winRatio >= _requiredSortedWinRatio) {
        highCount = mid;
        AutoPasLog(TRACE, "SortingThresholdBenchmark search {} cell direction={} n={}: sorted won {}/{} reps → high={}",
                   newton3, cellDirection, mid, outcome, _repetitions, highCount);
      } else {
        lowCount = mid + 1;
        AutoPasLog(TRACE,
                   "SortingThresholdBenchmark search {} cell direction={} n={}: sorted won only {}/{} reps → low={}",
                   newton3, cellDirection, mid, outcome, _repetitions, lowCount);
      }
    }
    AutoPasLog(DEBUG, "SortingThresholdBenchmark {} cell direction={} threshold={}", newton3, cellDirection, lowCount);
    // If the threshold does not converge we disable the threshold by setting it to infinity.
    if (lowCount == upperBound) {
      return std::numeric_limits<unsigned long>::max();
    }
    return lowCount;
  }
};
 
}  // namespace autopas