equality.h

Go to the documentation of this file.

#pragma once
 
#include <array>
#include <cmath>
#include <cstddef>
#include <utility>
 
#include "clustering/ndarray.h"
 
namespace clustering::math {

template <class T, std::size_t N, Layout LA, Layout LB>
bool arrayEqual(const NDArray<T, N, LA> &a, const NDArray<T, N, LB> &b) noexcept {
  std::array<std::size_t, N> shape{};
  std::size_t total = 1;
  for (std::size_t k = 0; k < N; ++k) {
    if (a.dim(k) != b.dim(k)) {
      return false;
    }
    shape[k] = a.dim(k);
    total *= shape[k];
  }
  if (total == 0) {
    return true;
  }
  std::array<std::size_t, N> idx{};
  for (std::size_t flat = 0; flat < total; ++flat) {
    const bool equal = [&]<std::size_t... Ks>(std::index_sequence<Ks...>) {
      return a(idx[Ks]...) == b(idx[Ks]...);
    }(std::make_index_sequence<N>{});
    if (!equal) {
      return false;
    }
    for (std::size_t k = N; k-- > 0;) {
      if (++idx[k] < shape[k]) {
        break;
      }
      idx[k] = 0;
    }
  }
  return true;
}

template <class T, std::size_t N, Layout LA, Layout LB>
bool allClose(const NDArray<T, N, LA> &a, const NDArray<T, N, LB> &b, T rtol = T(1e-5),
              T atol = T(1e-8)) noexcept {
  std::array<std::size_t, N> shape{};
  std::size_t total = 1;
  for (std::size_t k = 0; k < N; ++k) {
    if (a.dim(k) != b.dim(k)) {
      return false;
    }
    shape[k] = a.dim(k);
    total *= shape[k];
  }
  if (total == 0) {
    return true;
  }
  std::array<std::size_t, N> idx{};
  for (std::size_t flat = 0; flat < total; ++flat) {
    auto [va, vb] = [&]<std::size_t... Ks>(std::index_sequence<Ks...>) {
      return std::pair<T, T>{a(idx[Ks]...), b(idx[Ks]...)};
    }(std::make_index_sequence<N>{});
    const T diff = std::fabs(va - vb);
    const T tol = atol + (rtol * std::fabs(vb));
    if (!(diff <= tol)) {
      return false;
    }
    for (std::size_t k = N; k-- > 0;) {
      if (++idx[k] < shape[k]) {
        break;
      }
      idx[k] = 0;
    }
  }
  return true;
}
 
} // namespace clustering::math