ArrayMath.h

Go to the documentation of this file.

 
#pragma once
 
#include <array>
#include <cmath>
#include <numeric>
#include <sstream>
 
#include "autopas/utils/ContainerConcept.h"
#include "autopas/utils/Math.h"
 
namespace autopas::utils::ArrayMath {
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> add(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = a[d] + b[d];
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> sub(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = a[d] - b[d];
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> min(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = std::min<T>(a[d], b[d]);
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr bool less(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  bool result = true;
  for (std::size_t d = 0; d < SIZE; ++d) {
    result = result and (a[d] < b[d]);
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> max(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = std::max<T>(a[d], b[d]);
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> mul(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = a[d] * b[d];
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> div(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = a[d] / b[d];
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> addScalar(const std::array<T, SIZE> &a, T s) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = a[d] + s;
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> subScalar(const std::array<T, SIZE> &a, T s) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = a[d] - s;
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> mulScalar(const std::array<T, SIZE> &a, T s) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = a[d] * s;
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> divScalar(const std::array<T, SIZE> &a, T s) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = a[d] / s;
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> divScalar(T s, const std::array<T, SIZE> &a) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = s / a[d];
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr T dot(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  T result = 0.;
  for (std::size_t i = 0; i < SIZE; i++) {
    result += a[i] * b[i];
  }
  return result;
}
 
template <class T>
[[nodiscard]] constexpr std::array<T, 3> cross(const std::array<T, 3> &a, const std::array<T, 3> &b) {
  return {
      a[1] * b[2] - a[2] * b[1],
      a[2] * b[0] - a[0] * b[2],
      a[0] * b[1] - a[1] * b[0],
  };
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr T L2Norm(const std::array<T, SIZE> &a) {
  return std::sqrt(dot(a, a));
}
 
template <class T>
[[nodiscard]] constexpr typename T::value_type prod(const T &a) {
  return std::accumulate(a.cbegin(), a.cend(), static_cast<typename T::value_type>(1), std::multiplies<>());
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> abs(const std::array<T, SIZE> &a) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = std::abs(a[d]);
  }
  return result;
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> normalize(const std::array<T, SIZE> &a) {
  return mulScalar(a, static_cast<T>(1) / L2Norm(a));
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<T, SIZE> ceil(const std::array<T, SIZE> &a) {
  std::array<T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = std::ceil(a[d]);
  }
  return result;
}
 
template <typename target_T = int, typename float_T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<target_T, SIZE> floorAndCast(const std::array<float_T, SIZE> &a) {
  std::array<target_T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = static_cast<target_T>(std::floor(a[d]));
  }
  return result;
}
 
template <class target_T = int, class float_T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<target_T, SIZE> ceilAndCast(const std::array<float_T, SIZE> &a) {
  std::array<target_T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = static_cast<target_T>(std::ceil(a[d]));
  }
  return result;
}
 
template <ContainerType Container>
bool isInUlp(Container lhs, Container rhs, unsigned int ulpDistance = Math::MAX_ULP_DISTANCE) {
  return lhs.size() == rhs.size() && std::ranges::equal(lhs, rhs, [&ulpDistance](const auto &lhs, const auto &rhs) {
           return Math::isInUlp(lhs, rhs, ulpDistance);
         });
}
 
template <ContainerType Container>
bool isNearRel(Container lhs, Container rhs, double maxRelativeDifference = Math::EPSILON_RELATIVE_EQUALITY) {
  return lhs.size() == rhs.size() &&
         std::ranges::equal(lhs, rhs, [&maxRelativeDifference](const auto &lhs, const auto &rhs) {
           return Math::isNearRel(lhs, rhs, maxRelativeDifference);
         });
}
 
template <class T, std::size_t SIZE>
[[nodiscard]] bool isNearRel(const std::vector<std::array<T, SIZE>> &a, const std::vector<std::array<T, SIZE>> &b,
                             double maxRelativeDifference = Math::EPSILON_RELATIVE_EQUALITY) {
  const auto size = a.size();
  if (size != b.size()) {
    return false;
  }
  bool arraysAreNear = true;
  for (std::size_t i = 0; i < size; ++i) {
    arraysAreNear = arraysAreNear and utils::ArrayMath::isNearRel(a[i], b[i], maxRelativeDifference);
  }
  return arraysAreNear;
}
 
template <class T>
[[nodiscard]] bool isEqual(const std::vector<T> &a, const std::vector<T> &b) {
  if (a.size() != b.size()) {
    return false;
  }
  bool arraysAreEqual = true;
  for (std::size_t i = 0; i < a.size(); ++i) {
    arraysAreEqual = arraysAreEqual and (a[i] == b[i]);
  }
  return arraysAreEqual;
}
 
template <typename new_T, typename old_T, std::size_t SIZE>
[[nodiscard]] constexpr std::array<new_T, SIZE> staticCastArray(const std::array<old_T, SIZE> &a) {
  std::array<new_T, SIZE> result{};
  for (std::size_t d = 0; d < SIZE; ++d) {
    result[d] = static_cast<new_T>(a[d]);
  }
  return result;
}
 
// namespace for templated operators
inline namespace literals {
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> operator+(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  return add(a, b);
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> &operator+=(std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  for (std::size_t d = 0; d < SIZE; ++d) {
    a[d] += b[d];
  }
  return a;
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> operator-(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  return sub(a, b);
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> &operator-=(std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  for (std::size_t d = 0; d < SIZE; ++d) {
    a[d] -= b[d];
  }
  return a;
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> operator*(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  return mul(a, b);
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> &operator*=(std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  for (std::size_t d = 0; d < SIZE; ++d) {
    a[d] *= b[d];
  }
  return a;
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> operator/(const std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  return div(a, b);
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> operator/(const std::array<T, SIZE> &a, T b) {
  return divScalar(a, b);
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> operator/(T a, const std::array<T, SIZE> &b) {
  return divScalar(a, b);
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> &operator/=(std::array<T, SIZE> &a, const std::array<T, SIZE> &b) {
  for (std::size_t d = 0; d < SIZE; ++d) {
    a[d] /= b[d];
  }
  return a;
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> operator+(const std::array<T, SIZE> &a, T s) {
  return addScalar(a, s);
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> &operator+=(std::array<T, SIZE> &a, T s) {
  for (std::size_t d = 0; d < SIZE; ++d) {
    a[d] += s;
  }
  return a;
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> operator-(const std::array<T, SIZE> &a, T s) {
  return subScalar(a, s);
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> &operator-=(std::array<T, SIZE> &a, T s) {
  for (std::size_t d = 0; d < SIZE; ++d) {
    a[d] -= s;
  }
  return a;
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> operator*(const std::array<T, SIZE> &a, T s) {
  return mulScalar(a, s);
}
 
template <class T, std::size_t SIZE>
constexpr std::array<T, SIZE> &operator*=(std::array<T, SIZE> &a, T s) {
  for (std::size_t d = 0; d < SIZE; ++d) {
    a[d] *= s;
  }
  return a;
}
}  // namespace literals
 
template <class T, std::size_t SIZE>
double boxDistanceSquared(const std::array<T, SIZE> &aMin, const std::array<T, SIZE> &aMax,
                          const std::array<T, SIZE> &bMin, const std::array<T, SIZE> &bMax) {
  using namespace autopas::utils::ArrayMath::literals;
  using autopas::utils::ArrayMath::dot;
  using autopas::utils::ArrayMath::max;
 
  const auto aToB = max(std::array<T, SIZE>{}, aMin - bMax);
  const auto bToA = max(std::array<T, SIZE>{}, bMin - aMax);
 
  return dot(aToB, aToB) + dot(bToA, bToA);
}
}  // namespace autopas::utils::ArrayMath