LiveInfo.h

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <variant>
 
#include "autopas/containers/CellBlock3D.h"
#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/ParticleBinStructure.h"
#include "autopas/utils/Timer.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.
 
 private:
  static utils::ParticleBinStructure buildCellBinStructure(const std::array<double, 3> &domainSize,
                                                           const double interactionLength,
                                                           const std::array<double, 3> &boxMin,
                                                           const std::array<double, 3> &boxMax, double cutoff) {
    std::array<size_t, 3> cellsPerDim{};
    std::array<double, 3> cellLength{};
 
    for (int d = 0; d < 3; ++d) {
      // The number of cells is rounded down because the cells will be stretched to fit.
      // std::max to ensure there is at least one cell.
      cellsPerDim[d] = std::max(static_cast<size_t>(std::floor(domainSize[d] / interactionLength)), 1ul);
      cellLength[d] = domainSize[d] / static_cast<double>(cellsPerDim[d]);
    }
 
    return {cellsPerDim, cellLength, boxMin, boxMax, cutoff};
  }
 
  static utils::ParticleBinStructure buildParticleDependentBinStructure(const std::array<double, 3> &domainSize,
                                                                        const size_t numParticles,
                                                                        const std::array<double, 3> &boxMin,
                                                                        const std::array<double, 3> &boxMax,
                                                                        double cutoff) {
    using namespace autopas::utils::ArrayMath::literals;
 
    const auto targetNumberOfBins = std::max(std::ceil(static_cast<double>(numParticles) / 10.), 1.);
    const auto targetNumberOfBinsPerDim = std::cbrt(targetNumberOfBins);
    // This is probably not an integer, so we floor to get more than 10 particles per bin
    const auto numberOfBinsPerDim = static_cast<size_t>(std::floor(targetNumberOfBinsPerDim));
    const auto binDimensions = domainSize / static_cast<double>(numberOfBinsPerDim);
 
    return {numberOfBinsPerDim, binDimensions, boxMin, boxMax, cutoff};
  }
 
  static utils::ParticleBinStructure buildBlurredBinStructure(const std::array<double, 3> &domainSize,
                                                              const std::array<double, 3> &boxMin,
                                                              const std::array<double, 3> &boxMax, double cutoff) {
    using namespace autopas::utils::ArrayMath::literals;
 
    const auto binLength = domainSize / 3.;
 
    return {3, binLength, boxMin, boxMax, cutoff};
  }
 
 public:
  using InfoType = std::variant<bool, double, size_t, ContainerOption, TraversalOption, LoadEstimatorOption,
                                DataLayoutOption, Newton3Option>;
 
  LiveInfo() = default;
 
  LiveInfo(size_t numOwnedParticles, size_t numHaloParticles, double cutoff, double skin, double domainSizeX,
           double domainSizeY, double domainSizeZ, size_t particleSize, size_t threadCount, size_t rebuildFrequency,
           size_t numCells, size_t numEmptyCells, size_t minParticlesPerCell, size_t maxParticlesPerCell,
           size_t medianParticlesPerCell, size_t lowerQuartileParticlesPerCell, size_t upperQuartileParticlesPerCell,
           double meanParticlesPerCell, double particlesPerCellStdDev, double relativeParticlesPerCellStdDev,
           size_t estimatedNumNeighborInteractions, double particleDependentBinMaxDensity,
           double particleDependentBinDensityStdDev, size_t maxParticlesPerBlurredBin, size_t minParticlesPerBlurredBin,
           size_t medianParticlesPerBlurredBin, size_t lowerQuartileParticlesPerBlurredBin,
           size_t upperQuartileParticlesPerBlurredBin, double meanParticlesPerBlurredBin,
           double particlesPerBlurredBinStdDev, double relativeParticlesPerBlurredBinStdDev) {
    // Validate all inputs and log warnings if invalid
    if (cutoff <= 0) {
      AutoPasLog(WARN, "LiveInfo: cutoff must be greater than 0");
    }
    if (skin < 0) {
      AutoPasLog(WARN, "LiveInfo: skin must be non-negative");
    }
    if (domainSizeX <= 0 or domainSizeY <= 0 or domainSizeZ <= 0) {
      AutoPasLog(WARN, "LiveInfo: domain sizes must be greater than 0");
    }
    if (particleSize == 0) {
      AutoPasLog(WARN, "LiveInfo: particleSize must be greater than 0");
    }
    if (threadCount == 0) {
      AutoPasLog(WARN, "LiveInfo: threadCount must be greater than 0");
    }
    if (rebuildFrequency == 0) {
      AutoPasLog(WARN, "LiveInfo: rebuildFrequency must be greater than 0");
    }
    if (numCells == 0) {
      AutoPasLog(WARN, "LiveInfo: numCells must be greater than 0");
    }
    if (maxParticlesPerCell < minParticlesPerCell) {
      AutoPasLog(WARN, "LiveInfo: maxParticlesPerCell must be >= minParticlesPerCell");
    }
    if (meanParticlesPerCell < 0) {
      AutoPasLog(WARN, "LiveInfo: meanParticlesPerCell must be non-negative");
    }
    if (particlesPerCellStdDev < 0) {
      AutoPasLog(WARN, "LiveInfo: particlesPerCellStdDev must be non-negative");
    }
    if (relativeParticlesPerCellStdDev < 0) {
      AutoPasLog(WARN, "LiveInfo: relativeParticlesPerCellStdDev must be non-negative");
    }
    if (particleDependentBinMaxDensity < 0) {
      AutoPasLog(WARN, "LiveInfo: particleDependentBinMaxDensity must be non-negative");
    }
    if (particleDependentBinDensityStdDev < 0) {
      AutoPasLog(WARN, "LiveInfo: particleDependentBinDensityStdDev must be non-negative");
    }
    if (maxParticlesPerBlurredBin < minParticlesPerBlurredBin) {
      AutoPasLog(WARN, "LiveInfo: maxParticlesPerBlurredBin must be >= minParticlesPerBlurredBin");
    }
    if (meanParticlesPerBlurredBin < 0) {
      AutoPasLog(WARN, "LiveInfo: meanParticlesPerBlurredBin must be non-negative");
    }
    if (particlesPerBlurredBinStdDev < 0) {
      AutoPasLog(WARN, "LiveInfo: particlesPerBlurredBinStdDev must be non-negative");
    }
    if (relativeParticlesPerBlurredBinStdDev < 0) {
      AutoPasLog(WARN, "LiveInfo: relativeParticlesPerBlurredBinStdDev must be non-negative");
    }
 
    // Store validated parameters in infos
    infos["numOwnedParticles"] = numOwnedParticles;
    infos["numHaloParticles"] = numHaloParticles;
    infos["cutoff"] = cutoff;
    infos["skin"] = skin;
    infos["domainSizeX"] = domainSizeX;
    infos["domainSizeY"] = domainSizeY;
    infos["domainSizeZ"] = domainSizeZ;
    infos["particleSize"] = particleSize;
    infos["threadCount"] = threadCount;
    infos["rebuildFrequency"] = rebuildFrequency;
    infos["numCells"] = numCells;
    infos["numEmptyCells"] = numEmptyCells;
    infos["minParticlesPerCell"] = minParticlesPerCell;
    infos["maxParticlesPerCell"] = maxParticlesPerCell;
    infos["medianParticlesPerCell"] = medianParticlesPerCell;
    infos["lowerQuartileParticlesPerCell"] = lowerQuartileParticlesPerCell;
    infos["upperQuartileParticlesPerCell"] = upperQuartileParticlesPerCell;
    infos["meanParticlesPerCell"] = meanParticlesPerCell;
    infos["particlesPerCellStdDev"] = particlesPerCellStdDev;
    infos["relativeParticlesPerCellStdDev"] = relativeParticlesPerCellStdDev;
    infos["estimatedNumNeighborInteractions"] = estimatedNumNeighborInteractions;
    infos["particleDependentBinMaxDensity"] = particleDependentBinMaxDensity;
    infos["particleDependentBinDensityStdDev"] = particleDependentBinDensityStdDev;
    infos["maxParticlesPerBlurredBin"] = maxParticlesPerBlurredBin;
    infos["minParticlesPerBlurredBin"] = minParticlesPerBlurredBin;
    infos["medianParticlesPerBlurredBin"] = medianParticlesPerBlurredBin;
    infos["lowerQuartileParticlesPerBlurredBin"] = lowerQuartileParticlesPerBlurredBin;
    infos["upperQuartileParticlesPerBlurredBin"] = upperQuartileParticlesPerBlurredBin;
    infos["meanParticlesPerBlurredBin"] = meanParticlesPerBlurredBin;
    infos["particlesPerBlurredBinStdDev"] = particlesPerBlurredBinStdDev;
    infos["relativeParticlesPerBlurredBinStdDev"] = relativeParticlesPerBlurredBinStdDev;
  }
 
  template <class Particle_T>
  void gather(ContainerIterator<Particle_T, true, false> particleIter, size_t rebuildFrequency,
              size_t numOwnedParticles, std::array<double, 3> boxMin, std::array<double, 3> boxMax, double cutoff,
              double skin) {
    using namespace utils::ArrayMath::literals;
    using utils::ArrayMath::ceilAndCast;
 
    // Timer for debugging
    utils::Timer timerGatherLiveInfo;
    timerGatherLiveInfo.start();
 
    // Aliases and info of particle distribution independent information
    const auto interactionLength = cutoff + skin;
 
    infos["cutoff"] = cutoff;
    infos["skin"] = skin;
    infos["rebuildFrequency"] = rebuildFrequency;
    infos["particleSize"] = sizeof(Particle_T);
    infos["threadCount"] = static_cast<size_t>(autopas_get_max_threads());
 
    infos["numOwnedParticles"] = numOwnedParticles;
 
    const auto domainSize = boxMax - boxMin;
    infos["domainSizeX"] = domainSize[0];
    infos["domainSizeY"] = domainSize[1];
    infos["domainSizeZ"] = domainSize[2];
 
    // Build bin structures
    auto cellBinStruct = buildCellBinStructure(domainSize, interactionLength, boxMin, boxMax, cutoff);
    auto particleDependentBinStruct =
        buildParticleDependentBinStructure(domainSize, numOwnedParticles, boxMin, boxMax, cutoff);
    auto blurredBinStruct = buildBlurredBinStructure(domainSize, boxMin, boxMax, cutoff);
 
    infos["numCells"] = cellBinStruct.getNumberOfBins();
 
    // Count the number of owned particles per bin for each bin structure. Also include total count for halo particles.
    size_t numOwnedParticlesCount = 0;
    size_t numHaloParticlesCount = 0;
    for (; particleIter.isValid(); ++particleIter) {
      if (particleIter->isOwned()) {
        numOwnedParticlesCount++;
        const auto particlePos = particleIter->getR();
        if (utils::inBox(particlePos, boxMin, boxMax)) {
          cellBinStruct.countParticle(particlePos);
          particleDependentBinStruct.countParticle(particlePos);
          blurredBinStruct.countParticle(particlePos);
        }
      } else if (particleIter->isHalo()) {
        numHaloParticlesCount++;
      }
    }
 
    // Sanity Check
    if (numOwnedParticlesCount != numOwnedParticles) {
      AutoPasLog(ERROR,
                 "Number of owned particles tracked by AutoPas ({}) does not match number of owned particles "
                 "counted using the iterator ({}).",
                 numOwnedParticles, numOwnedParticlesCount);
    }
 
    infos["numOwnedParticles"] = numOwnedParticlesCount;
    infos["numHaloParticles"] = numHaloParticlesCount;
 
    // calculate statistics
    cellBinStruct.calculateStatistics();
    particleDependentBinStruct.calculateStatistics();
    blurredBinStruct.calculateStatistics();
 
    // write cellBinStruct statistics to live info
    infos["numEmptyCells"] = cellBinStruct.getNumEmptyBins();
    infos["maxParticlesPerCell"] = cellBinStruct.getMaxParticlesPerBin();
    infos["minParticlesPerCell"] = cellBinStruct.getMinParticlesPerBin();
    infos["medianParticlesPerCell"] = cellBinStruct.getMedianParticlesPerBin();
    infos["lowerQuartileParticlesPerCell"] = cellBinStruct.getLowerQuartileParticlesPerBin();
    infos["upperQuartileParticlesPerCell"] = cellBinStruct.getUpperQuartileParticlesPerBin();
    infos["meanParticlesPerCell"] = cellBinStruct.getMeanParticlesPerBin();
    infos["relativeParticlesPerCellStdDev"] = cellBinStruct.getRelStdDevParticlesPerBin();
    infos["particlesPerCellStdDev"] = cellBinStruct.getStdDevParticlesPerBin();
    infos["estimatedNumNeighborInteractions"] = cellBinStruct.getEstimatedNumberOfNeighborInteractions();
 
    // write particle dependent bin statistics to live info
    infos["particleDependentBinMaxDensity"] = particleDependentBinStruct.getMaxDensity();
    infos["particleDependentBinDensityStdDev"] = particleDependentBinStruct.getStdDevDensity();
 
    // write blurred bin statistics to live info
    infos["maxParticlesPerBlurredBin"] = blurredBinStruct.getMaxParticlesPerBin();
    infos["minParticlesPerBlurredBin"] = blurredBinStruct.getMinParticlesPerBin();
    infos["medianParticlesPerBlurredBin"] = blurredBinStruct.getMedianParticlesPerBin();
    infos["lowerQuartileParticlesPerBlurredBin"] = blurredBinStruct.getLowerQuartileParticlesPerBin();
    infos["upperQuartileParticlesPerBlurredBin"] = blurredBinStruct.getUpperQuartileParticlesPerBin();
    infos["meanParticlesPerBlurredBin"] = blurredBinStruct.getMeanParticlesPerBin();
    infos["relativeParticlesPerBlurredBinStdDev"] = blurredBinStruct.getRelStdDevParticlesPerBin();
    infos["particlesPerBlurredBinStdDev"] = blurredBinStruct.getStdDevParticlesPerBin();
 
    timerGatherLiveInfo.stop();
    AutoPasLog(DEBUG, "Gathering of LiveInfo took {} ns.", timerGatherLiveInfo.getTotalTime());
  }
 
  [[nodiscard]] const auto &get() const { return infos; }
 
  template <typename T>
  T get(const std::string &key) const {
    // Find the key in the map
    const auto it = infos.find(key);
 
    // If key is not found, log an error
    if (it == infos.end()) {
      AutoPasLog(ERROR, "Key '" + key + "' not found in infos map.");
    }
 
    // Use std::get<T> to extract the value of the desired type
    try {
      return std::get<T>(it->second);
    } catch (const std::bad_variant_access &e) {
      AutoPasLog(ERROR, "Type mismatch for key '" + key + "'. Requested type does not match the stored type.");
      return T{};
    }
  }
 
  [[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 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