<commit_before>// -----------------------------------------------------------------------------
//
// Copyright (C) The BioDynaMo Project.
// All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
//
// See the LICENSE file distributed with this work for details.
// See the NOTICE file distributed with this work for additional information
// regarding copyright ownership.
//
// -----------------------------------------------------------------------------

#ifndef CORE_GRID_H_
#define CORE_GRID_H_

#include <assert.h>
#include <omp.h>

#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <iostream>
#include <limits>
#include <memory>
#include <mutex>
#ifdef LINUX
#include <parallel/algorithm>
#endif  // LINUX
#include <utility>
#include <vector>

#include <morton/morton.h>

#include "core/container/fixed_size_vector.h"
#include "core/container/inline_vector.h"
#include "core/container/math_array.h"
#include "core/container/parallel_resize_vector.h"
#include "core/container/sim_object_vector.h"
#include "core/functor.h"
#include "core/param/param.h"
#include "core/resource_manager.h"
#include "core/util/log.h"
#include "core/util/spinlock.h"

namespace bdm {

/// FIFO data structure. Stores max N number of objects.
/// Adding more elements will overwrite the oldest ones.
template <typename T, uint64_t N>
class CircularBuffer {
 public:
  CircularBuffer() {
    for (uint64_t i = 0; i < N; i++) {
      data_[i] = T();
    }
  }

  void clear() {  // NOLINT
    position_ = 0;
    for (uint64_t i = 0; i < N; i++) {
      data_[i].clear();
    }
  }

  void push_back(const T& data) {  // NOLINT
    data_[position_] = data;
    Increment();
  }

  T& operator[](uint64_t idx) { return data_[(idx + position_) % N]; }

  const T& operator[](uint64_t idx) const {
    return data_[(idx + position_) % N];
  }

  /// Calling function `End` and afterwards `Increment` is equivalent to
  /// `push_back`, but avoids copying data
  T* End() { return &(data_[position_]); }

  void Increment() {
    position_++;
    position_ %= N;
  }

 private:
  T data_[N];
  uint64_t position_ = 0;
};

/// A class that represents Cartesian 3D grid
class Grid {
  // DisplacementOpCuda needs access to some Grid private members to reconstruct
  // the grid on GPU (same for DisplacementOpOpenCL)
  friend class DisplacementOpCuda;
  friend class DisplacementOpOpenCL;

 public:
  /// A single unit cube of the grid
  struct Box {
    Spinlock lock_;
    // std::atomic<bool> timestamp_;
    uint32_t timestamp_;
    /// start value of the linked list of simulation objects inside this box.
    /// Next element can be found at `successors_[start_]`
    SoHandle start_;
    /// length of the linked list (i.e. number of simulation objects)
    /// uint64_t, because sizeof(Box) = 16, for uint16_t and uint64_t
    uint16_t length_;

    Box() : timestamp_(0), start_(SoHandle()), length_(0) {}
    /// Copy Constructor required for boxes_.resize()
    /// Since box values will be overwritten afterwards it forwards to the
    /// default ctor
    Box(const Box& other) : Box() {}

    Box& operator=(const Box& other) {
      // start_ = other.start_.load(std::memory_order_relaxed);
      // length_ = other.length_.load(std::memory_order_relaxed);
      start_ = other.start_;
      length_ = other.length_;
      return *this;
    }

    bool IsEmpty(uint64_t grid_timestamp) const {
      return grid_timestamp != timestamp_;
    }

    /// @brief      Adds a simulation object to this box
    ///
    /// @param[in]  so       The object's identifier
    /// @param   AddObject   successors   The successors
    void AddObject(SoHandle soh, SimObjectVector<SoHandle>* successors,
                   Grid* grid) {
      std::lock_guard<Spinlock> lock_guard(lock_);

      if (timestamp_ != grid->timestamp_) {
        timestamp_ = grid->timestamp_;
        length_ = 1;
        start_ = soh;
      } else {
        length_++;
        (*successors)[soh] = start_;
        start_ = soh;
      }
    }

    /// An iterator that iterates over the cells in this box
    struct Iterator {
      Iterator(Grid* grid, const Box* box)
          : grid_(grid), current_value_(box->start_), countdown_(box->length_) {
        if (grid->timestamp_ != box->timestamp_) {
          countdown_ = 0;
        }
      }

      bool IsAtEnd() { return countdown_ <= 0; }

      Iterator& operator++() {
        countdown_--;
        if (countdown_ > 0) {
          current_value_ = grid_->successors_[current_value_];
        }
        return *this;
      }

      SoHandle operator*() const { return current_value_; }

      /// Pointer to the neighbor grid; for accessing the successor_ list
      Grid* grid_;
      /// The current simulation object to be considered
      SoHandle current_value_;
      /// The remain number of simulation objects to consider
      int countdown_ = 0;
    };

    Iterator begin() const {  // NOLINT
      return Iterator(Simulation::GetActive()->GetGrid(), this);
    }
  };

  /// An iterator that iterates over the boxes in this grid
  struct NeighborIterator {
    explicit NeighborIterator(
        const FixedSizeVector<const Box*, 27>& neighbor_boxes,
        uint64_t grid_timestamp)
        : neighbor_boxes_(neighbor_boxes),
          // start iterator from box 0
          box_iterator_(neighbor_boxes_[0]->begin()),
          grid_timestamp_(grid_timestamp) {
      // if first box is empty
      if (neighbor_boxes_[0]->IsEmpty(grid_timestamp)) {
        ForwardToNonEmptyBox(grid_timestamp);
      }
    }

    bool IsAtEnd() const { return is_end_; }

    SoHandle operator*() const { return *box_iterator_; }

    /// Version where empty neighbor boxes are allowed
    NeighborIterator& operator++() {
      ++box_iterator_;
      // if iterator of current box has come to an end, continue with next box
      if (box_iterator_.IsAtEnd()) {
        return ForwardToNonEmptyBox(grid_timestamp_);
      }
      return *this;
    }

   private:
    /// The 27 neighbor boxes that will be searched for simulation objects
    const FixedSizeVector<const Box*, 27>& neighbor_boxes_;
    /// The box that shall be considered to iterate over for finding simulation
    /// objects
    typename Box::Iterator box_iterator_;
    uint64_t grid_timestamp_;
    /// The id of the box to be considered (i.e. value between 0 - 26)
    uint16_t box_idx_ = 0;
    /// Flag to indicate that all the neighbor boxes have been searched through
    bool is_end_ = false;

    /// Forwards the iterator to the next non empty box and returns itself
    /// If there are no non empty boxes is_end_ is set to true
    NeighborIterator& ForwardToNonEmptyBox(uint64_t grid_timestamp) {
      // increment box id until non empty box has been found
      while (++box_idx_ < neighbor_boxes_.size()) {
        // box is empty or uninitialized (padding box) -> continue
        if (neighbor_boxes_[box_idx_]->IsEmpty(grid_timestamp)) {
          continue;
        }
        // a non-empty box has been found
        box_iterator_ = neighbor_boxes_[box_idx_]->begin();
        return *this;
      }
      // all remaining boxes have been empty; reached end
      is_end_ = true;
      return *this;
    }
  };

  /// Enum that determines the degree of adjacency in search neighbor boxes
  //  todo(ahmad): currently only kHigh is supported (hardcoded 26 several
  //  places)
  enum Adjacency {
    kLow,    /**< The closest 8  neighboring boxes */
    kMedium, /**< The closest 18  neighboring boxes */
    kHigh    /**< The closest 26  neighboring boxes */
  };

  Grid() {}

  Grid(Grid const&) = delete;
  void operator=(Grid const&) = delete;

  /// @brief      Initialize the grid with the given simulation objects
  /// @param[in]  adjacency    The adjacency (see #Adjacency)
  void Initialize(Adjacency adjacency = kHigh) {
    adjacency_ = adjacency;

    UpdateGrid();
    initialized_ = true;
  }

  virtual ~Grid() {}

  /// Clears the grid
  void ClearGrid() {
    box_length_ = 1;
    largest_object_size_ = 0;
    num_boxes_axis_ = {{0}};
    num_boxes_xy_ = 0;
    int32_t inf = std::numeric_limits<int32_t>::max();
    grid_dimensions_ = {inf, -inf, inf, -inf, inf, -inf};
    threshold_dimensions_ = {inf, -inf};
    successors_.clear();
    has_grown_ = false;
  }

  struct AssignToBoxesFunctor : public Functor<void, SimObject*, SoHandle> {
    AssignToBoxesFunctor(Grid* grid) : grid_(grid) {}

    void operator()(SimObject* sim_object, SoHandle soh) override {
      const auto& position = sim_object->GetPosition();
      auto idx = grid_->GetBoxIndex(position);
      auto box = grid_->GetBoxPointer(idx);
      box->AddObject(soh, &(grid_->successors_), grid_);
      sim_object->SetBoxIdx(idx);
    }

   private:
    Grid* grid_ = nullptr;
  };

  /// Updates the grid, as simulation objects may have moved, added or deleted
  void UpdateGrid() {
    auto* rm = Simulation::GetActive()->GetResourceManager();

    if (rm->GetNumSimObjects() != 0) {
      ClearGrid();
      timestamp_++;

      auto inf = Math::kInfinity;
      std::array<double, 6> tmp_dim = {{inf, -inf, inf, -inf, inf, -inf}};
      CalculateGridDimensions(&tmp_dim);
      RoundOffGridDimensions(tmp_dim);

      auto los = ceil(largest_object_size_);
      assert(los > 0 &&
             "The largest object size was found to be 0. Please check if your "
             "cells are correctly initialized.");
      box_length_ = los;

      for (int i = 0; i < 3; i++) {
        int dimension_length =
            grid_dimensions_[2 * i + 1] - grid_dimensions_[2 * i];
        int r = dimension_length % box_length_;
        // If the grid is not perfectly divisible along each dimension by the
        // resolution, extend the grid so that it is
        if (r != 0) {
          // std::abs for the case that box_length_ > dimension_length
          grid_dimensions_[2 * i + 1] += (box_length_ - r);
        } else {
          // Else extend the grid dimension with one row, because the outmost
          // object lies exactly on the border
          grid_dimensions_[2 * i + 1] += box_length_;
        }
      }

      // Pad the grid to avoid out of bounds check when search neighbors
      for (int i = 0; i < 3; i++) {
        grid_dimensions_[2 * i] -= box_length_;
        grid_dimensions_[2 * i + 1] += box_length_;
      }

      // Calculate how many boxes fit along each dimension
      for (int i = 0; i < 3; i++) {
        int dimension_length =
            grid_dimensions_[2 * i + 1] - grid_dimensions_[2 * i];
        assert((dimension_length % box_length_ == 0) &&
               "The grid dimensions are not a multiple of its box length");
        num_boxes_axis_[i] = dimension_length / box_length_;
      }

      num_boxes_xy_ = num_boxes_axis_[0] * num_boxes_axis_[1];
      auto total_num_boxes = num_boxes_xy_ * num_boxes_axis_[2];

      CheckGridGrowth();

      // resize boxes_
      if (boxes_.size() != total_num_boxes) {
        if (boxes_.capacity() < total_num_boxes) {
          boxes_.reserve(total_num_boxes * 2);
        }
        boxes_.resize(total_num_boxes);
      }

      successors_.reserve();

      // Assign simulation objects to boxes
      AssignToBoxesFunctor functor(this);
      rm->ApplyOnAllElementsParallelDynamic(1000, functor);
      auto* param = Simulation::GetActive()->GetParam();
      if (param->bound_space_) {
        int min = param->min_bound_;
        int max = param->max_bound_;
        threshold_dimensions_ = {min, max};
      }

      if (param->thread_safety_mechanism_ ==
          Param::ThreadSafetyMechanism::kAutomatic) {
        nb_mutex_builder_->Update();
      }
    } else {
      // There are no sim objects in this simulation
      auto* param = Simulation::GetActive()->GetParam();

      bool uninitialized = boxes_.size() == 0;
      if (uninitialized && param->bound_space_) {
        // Simulation has never had any simulation objects
        // Initialize grid dimensions with `Param::min_bound_` and
        // `Param::max_bound_`
        // This is required for the DiffusionGrid
        int min = param->min_bound_;
        int max = param->max_bound_;
        grid_dimensions_ = {min, max, min, max, min, max};
        threshold_dimensions_ = {min, max};
        has_grown_ = true;
      } else if (!uninitialized) {
        // all simulation objects have been removed in the last iteration
        // grid state remains the same, but we have to set has_grown_ to false
        // otherwise the DiffusionGrid will attempt to resize
        has_grown_ = false;
      } else {
        Log::Fatal(
            "Grid",
            "You tried to initialize an empty simulation without bound space. "
            "Therefore we cannot determine the size of the simulation space. "
            "Please add simulation objects, or set Param::bound_space_, "
            "Param::min_bound_, and Param::max_bound_.");
      }
    }
  }

  /// @brief      Calculates the squared euclidian distance between two points
  ///             in 3D
  ///
  /// @param[in]  pos1  Position of the first point
  /// @param[in]  pos2  Position of the second point
  ///
  /// @return     The distance between the two points
  ///
  inline double SquaredEuclideanDistance(const Double3& pos1,
                                         const Double3& pos2) const {
    const double dx = pos2[0] - pos1[0];
    const double dy = pos2[1] - pos1[1];
    const double dz = pos2[2] - pos1[2];
    return (dx * dx + dy * dy + dz * dz);
  }

  inline bool WithinSquaredEuclideanDistance(double squared_radius,
                                             const Double3& pos1,
                                             const Double3& pos2) const {
    const double dx = pos2[0] - pos1[0];
    const double dx2 = dx * dx;
    if (dx2 > squared_radius) {
      return false;
    }

    const double dy = pos2[1] - pos1[1];
    const double dy2_plus_dx2 = dy * dy + dx2;
    if (dy2_plus_dx2 > squared_radius) {
      return false;
    }

    const double dz = pos2[2] - pos1[2];
    const double distance = dz * dz + dy2_plus_dx2;
    return distance < squared_radius;
  }

  void UpdateBoxZOrder() {
    // iterate boxes in Z-order / morton order
    // TODO(lukas) this is a very quick attempt to test an idea
    // improve performance of this brute force solution
    zorder_sorted_boxes_.resize(boxes_.size());
#pragma omp parallel for collapse(3)
    for (uint32_t x = 0; x < num_boxes_axis_[0]; x++) {
      for (uint32_t y = 0; y < num_boxes_axis_[1]; y++) {
        for (uint32_t z = 0; z < num_boxes_axis_[2]; z++) {
          auto box_idx = GetBoxIndex(std::array<uint32_t, 3>{x, y, z});
          auto morton = libmorton::morton3D_64_encode(x, y, z);
          zorder_sorted_boxes_[box_idx] =
              std::pair<uint32_t, const Box*>{morton, &boxes_[box_idx]};
        }
      }
    }
#ifdef LINUX
    __gnu_parallel::sort(
        zorder_sorted_boxes_.begin(), zorder_sorted_boxes_.end(),
        [](const auto& lhs, const auto& rhs) { return lhs.first < rhs.first; });
#else
    std::sort(
        zorder_sorted_boxes_.begin(), zorder_sorted_boxes_.end(),
        [](const auto& lhs, const auto& rhs) { return lhs.first < rhs.first; });
#endif  // LINUX
  }

  /// This method iterates over all elements. Iteration is performed in
  /// Z-order of boxes. There is no particular order for elements inside a box.
  template <typename Lambda>
  void IterateZOrder(const Lambda& lambda) {
    UpdateBoxZOrder();
    for (uint64_t i = 0; i < zorder_sorted_boxes_.size(); i++) {<commit_msg>Fix compilation errors using Intel compilers<commit_after>// -----------------------------------------------------------------------------
//
// Copyright (C) The BioDynaMo Project.
// All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
//
// See the LICENSE file distributed with this work for details.
// See the NOTICE file distributed with this work for additional information
// regarding copyright ownership.
//
// -----------------------------------------------------------------------------

#ifndef CORE_GRID_H_
#define CORE_GRID_H_

#include <assert.h>
#include <omp.h>

#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <iostream>
#include <limits>
#include <memory>
#include <mutex>
#ifdef LINUX
#include <parallel/algorithm>
#endif  // LINUX
#include <utility>
#include <vector>

#include <morton/morton.h>

#include "core/container/fixed_size_vector.h"
#include "core/container/inline_vector.h"
#include "core/container/math_array.h"
#include "core/container/parallel_resize_vector.h"
#include "core/container/sim_object_vector.h"
#include "core/functor.h"
#include "core/param/param.h"
#include "core/resource_manager.h"
#include "core/util/log.h"
#include "core/util/spinlock.h"

namespace bdm {

/// FIFO data structure. Stores max N number of objects.
/// Adding more elements will overwrite the oldest ones.
template <typename T, uint64_t N>
class CircularBuffer {
 public:
  CircularBuffer() {
    for (uint64_t i = 0; i < N; i++) {
      data_[i] = T();
    }
  }

  void clear() {  // NOLINT
    position_ = 0;
    for (uint64_t i = 0; i < N; i++) {
      data_[i].clear();
    }
  }

  void push_back(const T& data) {  // NOLINT
    data_[position_] = data;
    Increment();
  }

  T& operator[](uint64_t idx) { return data_[(idx + position_) % N]; }

  const T& operator[](uint64_t idx) const {
    return data_[(idx + position_) % N];
  }

  /// Calling function `End` and afterwards `Increment` is equivalent to
  /// `push_back`, but avoids copying data
  T* End() { return &(data_[position_]); }

  void Increment() {
    position_++;
    position_ %= N;
  }

 private:
  T data_[N];
  uint64_t position_ = 0;
};

/// A class that represents Cartesian 3D grid
class Grid {
  // DisplacementOpCuda needs access to some Grid private members to reconstruct
  // the grid on GPU (same for DisplacementOpOpenCL)
  friend class DisplacementOpCuda;
  friend class DisplacementOpOpenCL;

 public:
  /// A single unit cube of the grid
  struct Box {
    Spinlock lock_;
    // std::atomic<bool> timestamp_;
    uint32_t timestamp_;
    /// start value of the linked list of simulation objects inside this box.
    /// Next element can be found at `successors_[start_]`
    SoHandle start_;
    /// length of the linked list (i.e. number of simulation objects)
    /// uint64_t, because sizeof(Box) = 16, for uint16_t and uint64_t
    uint16_t length_;

    Box() : timestamp_(0), start_(SoHandle()), length_(0) {}
    /// Copy Constructor required for boxes_.resize()
    /// Since box values will be overwritten afterwards it forwards to the
    /// default ctor
    Box(const Box& other) : Box() {}

    Box& operator=(const Box& other) {
      // start_ = other.start_.load(std::memory_order_relaxed);
      // length_ = other.length_.load(std::memory_order_relaxed);
      start_ = other.start_;
      length_ = other.length_;
      return *this;
    }

    bool IsEmpty(uint64_t grid_timestamp) const {
      return grid_timestamp != timestamp_;
    }

    /// @brief      Adds a simulation object to this box
    ///
    /// @param[in]  so       The object's identifier
    /// @param   AddObject   successors   The successors
    void AddObject(SoHandle soh, SimObjectVector<SoHandle>* successors,
                   Grid* grid) {
      std::lock_guard<Spinlock> lock_guard(lock_);

      if (timestamp_ != grid->timestamp_) {
        timestamp_ = grid->timestamp_;
        length_ = 1;
        start_ = soh;
      } else {
        length_++;
        (*successors)[soh] = start_;
        start_ = soh;
      }
    }

    /// An iterator that iterates over the cells in this box
    struct Iterator {
      Iterator(Grid* grid, const Box* box)
          : grid_(grid), current_value_(box->start_), countdown_(box->length_) {
        if (grid->timestamp_ != box->timestamp_) {
          countdown_ = 0;
        }
      }

      bool IsAtEnd() { return countdown_ <= 0; }

      Iterator& operator++() {
        countdown_--;
        if (countdown_ > 0) {
          current_value_ = grid_->successors_[current_value_];
        }
        return *this;
      }

      SoHandle operator*() const { return current_value_; }

      /// Pointer to the neighbor grid; for accessing the successor_ list
      Grid* grid_;
      /// The current simulation object to be considered
      SoHandle current_value_;
      /// The remain number of simulation objects to consider
      int countdown_ = 0;
    };

    Iterator begin() const {  // NOLINT
      return Iterator(Simulation::GetActive()->GetGrid(), this);
    }
  };

  /// An iterator that iterates over the boxes in this grid
  struct NeighborIterator {
    explicit NeighborIterator(
        const FixedSizeVector<const Box*, 27>& neighbor_boxes,
        uint64_t grid_timestamp)
        : neighbor_boxes_(neighbor_boxes),
          // start iterator from box 0
          box_iterator_(neighbor_boxes_[0]->begin()),
          grid_timestamp_(grid_timestamp) {
      // if first box is empty
      if (neighbor_boxes_[0]->IsEmpty(grid_timestamp)) {
        ForwardToNonEmptyBox(grid_timestamp);
      }
    }

    bool IsAtEnd() const { return is_end_; }

    SoHandle operator*() const { return *box_iterator_; }

    /// Version where empty neighbor boxes are allowed
    NeighborIterator& operator++() {
      ++box_iterator_;
      // if iterator of current box has come to an end, continue with next box
      if (box_iterator_.IsAtEnd()) {
        return ForwardToNonEmptyBox(grid_timestamp_);
      }
      return *this;
    }

   private:
    /// The 27 neighbor boxes that will be searched for simulation objects
    const FixedSizeVector<const Box*, 27>& neighbor_boxes_;
    /// The box that shall be considered to iterate over for finding simulation
    /// objects
    typename Box::Iterator box_iterator_;
    uint64_t grid_timestamp_;
    /// The id of the box to be considered (i.e. value between 0 - 26)
    uint16_t box_idx_ = 0;
    /// Flag to indicate that all the neighbor boxes have been searched through
    bool is_end_ = false;

    /// Forwards the iterator to the next non empty box and returns itself
    /// If there are no non empty boxes is_end_ is set to true
    NeighborIterator& ForwardToNonEmptyBox(uint64_t grid_timestamp) {
      // increment box id until non empty box has been found
      while (++box_idx_ < neighbor_boxes_.size()) {
        // box is empty or uninitialized (padding box) -> continue
        if (neighbor_boxes_[box_idx_]->IsEmpty(grid_timestamp)) {
          continue;
        }
        // a non-empty box has been found
        box_iterator_ = neighbor_boxes_[box_idx_]->begin();
        return *this;
      }
      // all remaining boxes have been empty; reached end
      is_end_ = true;
      return *this;
    }
  };

  /// Enum that determines the degree of adjacency in search neighbor boxes
  //  todo(ahmad): currently only kHigh is supported (hardcoded 26 several
  //  places)
  enum Adjacency {
    kLow,    /**< The closest 8  neighboring boxes */
    kMedium, /**< The closest 18  neighboring boxes */
    kHigh    /**< The closest 26  neighboring boxes */
  };

  Grid() {}

  Grid(Grid const&) = delete;
  void operator=(Grid const&) = delete;

  /// @brief      Initialize the grid with the given simulation objects
  /// @param[in]  adjacency    The adjacency (see #Adjacency)
  void Initialize(Adjacency adjacency = kHigh) {
    adjacency_ = adjacency;

    UpdateGrid();
    initialized_ = true;
  }

  virtual ~Grid() {}

  /// Clears the grid
  void ClearGrid() {
    box_length_ = 1;
    largest_object_size_ = 0;
    num_boxes_axis_ = {{0}};
    num_boxes_xy_ = 0;
    int32_t inf = std::numeric_limits<int32_t>::max();
    grid_dimensions_ = {inf, -inf, inf, -inf, inf, -inf};
    threshold_dimensions_ = {inf, -inf};
    successors_.clear();
    has_grown_ = false;
  }

  struct AssignToBoxesFunctor : public Functor<void, SimObject*, SoHandle> {
    AssignToBoxesFunctor(Grid* grid) : grid_(grid) {}

    void operator()(SimObject* sim_object, SoHandle soh) override {
      const auto& position = sim_object->GetPosition();
      auto idx = grid_->GetBoxIndex(position);
      auto box = grid_->GetBoxPointer(idx);
      box->AddObject(soh, &(grid_->successors_), grid_);
      sim_object->SetBoxIdx(idx);
    }

   private:
    Grid* grid_ = nullptr;
  };

  /// Updates the grid, as simulation objects may have moved, added or deleted
  void UpdateGrid() {
    auto* rm = Simulation::GetActive()->GetResourceManager();

    if (rm->GetNumSimObjects() != 0) {
      ClearGrid();
      timestamp_++;

      auto inf = Math::kInfinity;
      std::array<double, 6> tmp_dim = {{inf, -inf, inf, -inf, inf, -inf}};
      CalculateGridDimensions(&tmp_dim);
      RoundOffGridDimensions(tmp_dim);

      auto los = ceil(largest_object_size_);
      assert(los > 0 &&
             "The largest object size was found to be 0. Please check if your "
             "cells are correctly initialized.");
      box_length_ = los;

      for (int i = 0; i < 3; i++) {
        int dimension_length =
            grid_dimensions_[2 * i + 1] - grid_dimensions_[2 * i];
        int r = dimension_length % box_length_;
        // If the grid is not perfectly divisible along each dimension by the
        // resolution, extend the grid so that it is
        if (r != 0) {
          // std::abs for the case that box_length_ > dimension_length
          grid_dimensions_[2 * i + 1] += (box_length_ - r);
        } else {
          // Else extend the grid dimension with one row, because the outmost
          // object lies exactly on the border
          grid_dimensions_[2 * i + 1] += box_length_;
        }
      }

      // Pad the grid to avoid out of bounds check when search neighbors
      for (int i = 0; i < 3; i++) {
        grid_dimensions_[2 * i] -= box_length_;
        grid_dimensions_[2 * i + 1] += box_length_;
      }

      // Calculate how many boxes fit along each dimension
      for (int i = 0; i < 3; i++) {
        int dimension_length =
            grid_dimensions_[2 * i + 1] - grid_dimensions_[2 * i];
        assert((dimension_length % box_length_ == 0) &&
               "The grid dimensions are not a multiple of its box length");
        num_boxes_axis_[i] = dimension_length / box_length_;
      }

      num_boxes_xy_ = num_boxes_axis_[0] * num_boxes_axis_[1];
      auto total_num_boxes = num_boxes_xy_ * num_boxes_axis_[2];

      CheckGridGrowth();

      // resize boxes_
      if (boxes_.size() != total_num_boxes) {
        if (boxes_.capacity() < total_num_boxes) {
          boxes_.reserve(total_num_boxes * 2);
        }
        boxes_.resize(total_num_boxes);
      }

      successors_.reserve();

      // Assign simulation objects to boxes
      AssignToBoxesFunctor functor(this);
      rm->ApplyOnAllElementsParallelDynamic(1000, functor);
      auto* param = Simulation::GetActive()->GetParam();
      if (param->bound_space_) {
        int min = param->min_bound_;
        int max = param->max_bound_;
        threshold_dimensions_ = {min, max};
      }

      if (param->thread_safety_mechanism_ ==
          Param::ThreadSafetyMechanism::kAutomatic) {
        nb_mutex_builder_->Update();
      }
    } else {
      // There are no sim objects in this simulation
      auto* param = Simulation::GetActive()->GetParam();

      bool uninitialized = boxes_.size() == 0;
      if (uninitialized && param->bound_space_) {
        // Simulation has never had any simulation objects
        // Initialize grid dimensions with `Param::min_bound_` and
        // `Param::max_bound_`
        // This is required for the DiffusionGrid
        int min = param->min_bound_;
        int max = param->max_bound_;
        grid_dimensions_ = {min, max, min, max, min, max};
        threshold_dimensions_ = {min, max};
        has_grown_ = true;
      } else if (!uninitialized) {
        // all simulation objects have been removed in the last iteration
        // grid state remains the same, but we have to set has_grown_ to false
        // otherwise the DiffusionGrid will attempt to resize
        has_grown_ = false;
      } else {
        Log::Fatal(
            "Grid",
            "You tried to initialize an empty simulation without bound space. "
            "Therefore we cannot determine the size of the simulation space. "
            "Please add simulation objects, or set Param::bound_space_, "
            "Param::min_bound_, and Param::max_bound_.");
      }
    }
  }

  /// @brief      Calculates the squared euclidian distance between two points
  ///             in 3D
  ///
  /// @param[in]  pos1  Position of the first point
  /// @param[in]  pos2  Position of the second point
  ///
  /// @return     The distance between the two points
  ///
  inline double SquaredEuclideanDistance(const Double3& pos1,
                                         const Double3& pos2) const {
    const double dx = pos2[0] - pos1[0];
    const double dy = pos2[1] - pos1[1];
    const double dz = pos2[2] - pos1[2];
    return (dx * dx + dy * dy + dz * dz);
  }

  inline bool WithinSquaredEuclideanDistance(double squared_radius,
                                             const Double3& pos1,
                                             const Double3& pos2) const {
    const double dx = pos2[0] - pos1[0];
    const double dx2 = dx * dx;
    if (dx2 > squared_radius) {
      return false;
    }

    const double dy = pos2[1] - pos1[1];
    const double dy2_plus_dx2 = dy * dy + dx2;
    if (dy2_plus_dx2 > squared_radius) {
      return false;
    }

    const double dz = pos2[2] - pos1[2];
    const double distance = dz * dz + dy2_plus_dx2;
    return distance < squared_radius;
  }

  void UpdateBoxZOrder() {
    // iterate boxes in Z-order / morton order
    // TODO(lukas) this is a very quick attempt to test an idea
    // improve performance of this brute force solution
    zorder_sorted_boxes_.resize(boxes_.size());
    const uint32_t nx = num_boxes_axis_[0];
    const uint32_t ny = num_boxes_axis_[1];
    const uint32_t nz = num_boxes_axis_[2];
#pragma omp parallel for collapse(3)
    for (uint32_t x = 0; x < nx; x++) {
      for (uint32_t y = 0; y < ny; y++) {
        for (uint32_t z = 0; z < nz; z++) {
          auto box_idx = GetBoxIndex(std::array<uint32_t, 3>{x, y, z});
          auto morton = libmorton::morton3D_64_encode(x, y, z);
          zorder_sorted_boxes_[box_idx] =
              std::pair<uint32_t, const Box*>{morton, &boxes_[box_idx]};
        }
      }
    }
#ifdef LINUX
    __gnu_parallel::sort(
        zorder_sorted_boxes_.begin(), zorder_sorted_boxes_.end(),
        [](const auto& lhs, const auto& rhs) { return lhs.first < rhs.first; });
#else
    std::sort(
        zorder_sorted_boxes_.begin(), zorder_sorted_boxes_.end(),
        [](const auto& lhs, const auto& rhs) { return lhs.first < rhs.first; });
#endif  // LINUX
  }

  /// This method iterates over all elements. Iteration is performed in
  /// Z-order of boxes. There is no particular order for elements inside a box.
  template <typename Lambda>
  void IterateZOrder(const Lambda& lambda) {
    UpdateBoxZOrder();
    for (uint64_t i = 0; i < zorder_sorted_boxes_.size(); i++) {
{'fraction_non_alphanumeric': 0.11909520637035627, 'fraction_numerical': 0.011484407220583104, 'mean_word_length': 2.4236993591832303, 'pattern_counts': {'":': 0, '<': 91, '<?xml version=': 0, '>': 129, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1405557', 'n_tokens_mistral': 9729, 'n_tokens_neox': 9299, 'n_words': 2929}