LiveInfo.h

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <variant>
 
#include "autopas/containers/ParticleContainerInterface.h"
#include "autopas/options/ContainerOption.h"
#include "autopas/options/DataLayoutOption.h"
#include "autopas/options/LoadEstimatorOption.h"
#include "autopas/options/Newton3Option.h"
#include "autopas/options/TraversalOption.h"
#include "autopas/utils/ArrayMath.h"
#include "autopas/utils/SimilarityFunctions.h"
#include "autopas/utils/WrapOpenMP.h"
 
namespace autopas {
 
class LiveInfo {
  // The class currently needs to be defined in a header only since iteration over a particle container requires to
  // know the Particle type. Actually, no particle specific information can be used here, so only a pointer to
  // ParticleBase would suffice, but iteration doesn't work that way at the moment.
 
 public:
  using InfoType = std::variant<bool, double, size_t, ContainerOption, TraversalOption, LoadEstimatorOption,
                                DataLayoutOption, Newton3Option>;
 
  template <class Particle_T, class PairwiseFunctor>
  void gather(const autopas::ParticleContainerInterface<Particle_T> &container, const PairwiseFunctor &functor,
              unsigned int rebuildFrequency) {
    using namespace autopas::utils::ArrayMath::literals;
    using autopas::utils::ArrayMath::ceilToInt;
    using autopas::utils::ArrayMath::floorToInt;
 
    // Some aliases for quicker access
    const auto &boxMin = container.getBoxMin();
    const auto &boxMax = container.getBoxMax();
    const auto cutoff = container.getCutoff();
    const auto cutoffInv = 1.0 / cutoff;
    const auto numParticles = container.getNumberOfParticles();
 
    infos["numParticles"] = numParticles;
    infos["cutoff"] = cutoff;
    infos["skin"] = container.getVerletSkin();
    infos["rebuildFrequency"] = static_cast<size_t>(rebuildFrequency);
    const auto domainSize = boxMax - boxMin;
 
    infos["domainSizeX"] = domainSize[0];
    infos["domainSizeY"] = domainSize[1];
    infos["domainSizeZ"] = domainSize[2];
 
    infos["particleSize"] = sizeof(Particle_T);
 
    // Calculate number of cells for a linked cells container, assuming cell size factor == 1
    const auto cellsPerDim = ceilToInt(domainSize * cutoffInv);
    const auto numCells = static_cast<size_t>(cellsPerDim[0] * cellsPerDim[1] * cellsPerDim[2]);
    infos["numCells"] = numCells;
 
    // Count how many particles are in each cell via bin counting
    std::vector<size_t> particleBins;
    // +1 because we count all halo particles in the last bin
    particleBins.resize(numCells + 1);
    // Blurred cells divide the domain into 3x3x3 equivalent boxes.
    std::vector<size_t> particleBinsBlurred;
    particleBinsBlurred.resize(27);
    const auto blurredCellDimsReciproc = std::array<double, 3>{3.0, 3.0, 3.0} / domainSize;
    for (const Particle_T &particle : container) {
      if (utils::inBox(particle.getR(), boxMin, boxMax)) {
        // find the actual cell
        const auto offsetIntoBox = particle.getR() - boxMin;
        const auto cell = floorToInt(offsetIntoBox * cutoffInv);
        const auto binIndex = (cell[2] * cellsPerDim[1] + cell[1]) * cellsPerDim[0] + cell[0];
        particleBins[binIndex]++;
 
        // find the blurred cell
        const auto cellBlurred = floorToInt(offsetIntoBox * blurredCellDimsReciproc);
        const auto binIndexBlurred = (cellBlurred[2] * 3 + cellBlurred[1]) * 3 + cellBlurred[0];
        particleBinsBlurred[binIndexBlurred]++;
      } else {
        // found a halo particle
        particleBins.back()++;
      }
    }
 
    infos["numHaloParticles"] = particleBins.back();
 
    // calculate statistics about particle distributions per cell
    const auto avgParticlesPerCell = numParticles == 0
                                         ? 0.
                                         : static_cast<double>(container.getNumberOfParticles() - particleBins.back()) /
                                               static_cast<double>(numCells);
    const auto avgParticlesPerBlurredCell =
        numParticles == 0 ? 0.
                          : static_cast<double>(container.getNumberOfParticles() - particleBins.back()) /
                                static_cast<double>(particleBinsBlurred.size());
 
    const auto [estimatedNumNeighborInteractions, maxDiff, sumStddev, numEmptyCells, maxParticlesPerCell,
                minParticlesPerCell] = [&]() {
      double estimatedNumNeighborInteractionsLambda = 0.;
      double maxDiffLambda = 0.;
      double sumStddevLambda = 0.;
      size_t numEmptyCellsLambda = 0;
      size_t maxParticlesPerCellLambda = std::numeric_limits<size_t>::min();
      size_t minParticlesPerCellLambda = std::numeric_limits<size_t>::max();
      // go over all bins and calculate statistics
      for (size_t i = 0; i < particleBins.size() - 1; i++) {
        const auto particlesInBin = particleBins[i];
        if (particlesInBin == 0) {
          ++numEmptyCellsLambda;
        } else {
          const auto estimatedNumParticlesInVicinity = particlesInBin * (particlesInBin * 27 - 1);
          // ratio of sphere volume to volume of 3x3x3 cubes where edge size of cube is same as the radius of sphere
          constexpr double probabilityParticleWithinCutoff = 0.155;
          estimatedNumNeighborInteractionsLambda += estimatedNumParticlesInVicinity * probabilityParticleWithinCutoff;
        }
        maxParticlesPerCellLambda = std::max(particlesInBin, maxParticlesPerCellLambda);
        minParticlesPerCellLambda = std::min(particlesInBin, minParticlesPerCellLambda);
        const auto diffFromAvg = avgParticlesPerCell - static_cast<int>(particlesInBin);
        maxDiffLambda = std::max(diffFromAvg, maxDiffLambda);
        sumStddevLambda += diffFromAvg * diffFromAvg;
      }
      estimatedNumNeighborInteractionsLambda /= 2;
 
      return std::tuple{estimatedNumNeighborInteractionsLambda,
                        maxDiffLambda,
                        sumStddevLambda,
                        numEmptyCellsLambda,
                        maxParticlesPerCellLambda,
                        minParticlesPerCellLambda};
    }();
 
    infos["numEmptyCells"] = numEmptyCells;
    infos["maxParticlesPerCell"] = maxParticlesPerCell;
    infos["minParticlesPerCell"] = minParticlesPerCell;
    infos["particlesPerCellStdDev"] =
        numParticles == 0 ? 0.
                          : std::sqrt(sumStddev) / static_cast<double>(particleBins.size() - 1) / avgParticlesPerCell;
    infos["avgParticlesPerCell"] = avgParticlesPerCell;
 
    const auto [homogeneity, maxDensity] = autopas::utils::calculateHomogeneityAndMaxDensity(container);
 
    infos["homogeneity"] = homogeneity;
    infos["maxDensity"] = maxDensity;
 
    infos["estimatedNumNeighborInteractions"] = static_cast<unsigned long>(estimatedNumNeighborInteractions);
 
    const double sumStddevBlurred = [&]() {
      double res = 0.;
      for (auto numParticlesInBin : particleBinsBlurred) {
        auto diff = avgParticlesPerBlurredCell - static_cast<int>(numParticlesInBin);
        res += diff * diff;
      }
      return res;
    }();
    infos["particlesPerBlurredCellStdDev"] = numParticles == 0 ? 0.
                                                               : std::sqrt(sumStddevBlurred) /
                                                                     static_cast<double>(particleBinsBlurred.size()) /
                                                                     avgParticlesPerBlurredCell;
 
    infos["threadCount"] = static_cast<size_t>(autopas::autopas_get_max_threads());
 
    constexpr size_t particleSizeNeededByFunctor = calculateParticleSizeNeededByFunctor<Particle_T, PairwiseFunctor>(
        std::make_index_sequence<PairwiseFunctor::getNeededAttr().size()>());
    infos["particleSizeNeededByFunctor"] = particleSizeNeededByFunctor;
  }
 
  [[nodiscard]] const auto &get() const { return infos; }
 
  [[nodiscard]] std::string toString() const {
    auto typeToString = [](auto type) {
      if constexpr (std::is_same_v<decltype(type), bool> or std::is_same_v<decltype(type), double> or
                    std::is_same_v<decltype(type), size_t>) {
        return std::to_string(type);
      } else if constexpr (std::is_base_of_v<Option<decltype(type)>, decltype(type)>) {
        return type.to_string();
      }
      return std::string{"fail"};
    };
    auto pairToString = [&](const auto &pair) { return pair.first + "=" + std::visit(typeToString, pair.second); };
    std::string res{"Live Info: "};
    if (not infos.empty()) {
      // We initialize the accumulation with the first element. Hence, the accumulation starts at next(begin).
      res += std::accumulate(std::next(infos.begin()), infos.end(), pairToString(*infos.begin()),
                             [&](std::string s, const auto &elem) { return std::move(s) + " " + pairToString(elem); });
    }
    return res;
  }
 
  friend std::ostream &operator<<(std::ostream &out, const LiveInfo &info) {
    out << info.infos.size() << ' ';
    for (const auto &[name, val] : info.infos) {
      // val.index here is the index of this value's type in the LiveInfo::InfoType variant type.
      // This is needed for determining the type when reading this file via readIndex.
      out << name << ' ' << val.index() << ' ';
      std::visit([&](const auto &v) { out << v << ' '; }, val);
    }
    return out;
  }
 
  friend std::istream &operator>>(std::istream &in, LiveInfo &info) {
    size_t numElements{0};
    in >> numElements;
    for (size_t i = 0; i < numElements; i++) {
      std::string name;
      in >> name;
      size_t idx{0};
      in >> idx;
      const auto val =
          readIndex<LiveInfo::InfoType>(in, idx, std::make_index_sequence<std::variant_size_v<LiveInfo::InfoType>>());
      info.infos[name] = val;
    }
    return in;
  }
 
  [[nodiscard]] std::pair<std::string, std::string> getCSVLine() const {
    // match all words (anything that is neither a ' ' or '='), that are followed by a '=',
    // ignoring the 'Live Info: ' prefix
    const auto keyRegex = std::regex("([^= ]+)=[^ ]*");
    // match all words that are preceded by a '='
    const auto valueRegex = std::regex("=([^ ]+)");
 
    auto searchString = toString();
    // remove leading Live Info:
    searchString = searchString.substr(std::string("Live Info: ").size());
 
    std::sregex_iterator keyIter(searchString.begin(), searchString.end(), keyRegex);
    std::sregex_iterator valueIter(searchString.begin(), searchString.end(), valueRegex);
    std::sregex_iterator end;
 
    std::stringstream header;
    std::stringstream line;
 
    while (keyIter != end) {
      // first submatch is the match of the capture group
      header << keyIter->str(1) << ",";
      ++keyIter;
    }
    while (valueIter != end) {
      // first submatch is the match of the capture group
      line << valueIter->str(1) << ",";
      ++valueIter;
    }
 
    auto headerStr = header.str();
    auto lineStr = line.str();
    // drop trailing ','
    headerStr.pop_back();
    lineStr.pop_back();
 
    return std::make_pair(headerStr, lineStr);
  }
 
 private:
  template <class Particle_T, class PairwiseFunctor, size_t... Idx>
  constexpr static auto calculateParticleSizeNeededByFunctor(std::index_sequence<Idx...>) {
    return (0 + ... +
            sizeof(typename std::tuple_element<PairwiseFunctor::getNeededAttr()[Idx],
                                               typename Particle_T::SoAArraysType>::type::value_type));
  }
 
  template <class Variant, class Type, size_t Idx>
  static void readIndexHelper(std::istream &in, size_t idx, Variant &var) {
    if (Idx == idx) {
      Type val;
      in >> val;
      var = val;
    }
  }
 
  template <class Variant, size_t... Idx>
  static Variant readIndex(std::istream &in, size_t idx, std::index_sequence<Idx...>) {
    Variant var;
    (readIndexHelper<Variant, std::variant_alternative_t<Idx, Variant>, Idx>(in, idx, var), ...);
    return var;
  }
 
  std::map<std::string, InfoType> infos;
};
 
}  // namespace autopas