tpl_ca_neighborhood.H

Go to the documentation of this file.

/*
                          Aleph_w
 
  Data structures & Algorithms
  version 2.0.0b
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/
 
#ifndef TPL_CA_NEIGHBORHOOD_H
#define TPL_CA_NEIGHBORHOOD_H
 
#include <array>
#include <cstddef>
#include <span>
#include <utility>
 
#include <ah-errors.H>
#include <ca-traits.H>
#include <tpl_ca_concepts.H>
 
namespace Aleph {
namespace CA {
 
namespace ca_neighborhood_detail {
constexpr std::size_t static_pow(std::size_t base, std::size_t exp) noexcept
{
  std::size_t r = 1;
  for (std::size_t i = 0; i < exp; ++i)
    r *= base;
  return r;
}
 
template <std::size_t N, std::size_t R>
constexpr std::size_t von_neumann_count() noexcept
{
  const std::size_t total = static_pow(2 * R + 1, N);
  const ca_index_t Ri = static_cast<ca_index_t>(R);
 
  // Iterate over the Moore box and count points whose L1 norm lies
  // in [1, R]. Compile-time only; tiny N and R in practice.
  std::size_t cnt = 0;
  // current offset, initialised to (-R, -R, ..., -R).
  std::array<std::size_t, N> cur_buf{};
  for (std::size_t d = 0; d < N; ++d)
    cur_buf[d] = 0;  // we'll bias by -R when reading
 
  for (std::size_t i = 0; i < total; ++i)
    {
      // Convert cur_buf (in [0, 2R]) into signed offset in [-R, R].
      std::size_t l1 = 0;
      for (std::size_t d = 0; d < N; ++d)
        {
          const ca_index_t v = static_cast<ca_index_t>(cur_buf[d]) - Ri;
          l1 += static_cast<std::size_t>(v < 0 ? -v : v);
        }
      if (l1 != 0 and l1 <= R)
        ++cnt;
 
      // Increment cur_buf as a (2R+1)-base counter (last axis fastest).
      for (std::size_t d = N; d-- > 0;)
        {
          if (cur_buf[d] + 1 < 2 * R + 1)
            {
              ++cur_buf[d];
              break;
            }
          cur_buf[d] = 0;
        }
    }
  return cnt;
}
 
template <std::size_t N, std::size_t R>
constexpr auto compute_moore_offsets() noexcept
{
  constexpr std::size_t total = static_pow(2 * R + 1, N);
  constexpr std::size_t out_n = total - 1;
  std::array<Offset_Vec<N>, out_n> result{};
 
  Offset_Vec<N> cur{};
  for (std::size_t d = 0; d < N; ++d)
    cur[d] = -static_cast<ca_index_t>(R);
 
  std::size_t idx = 0;
  for (std::size_t i = 0; i < total; ++i)
    {
      bool zero = true;
      for (std::size_t d = 0; d < N; ++d)
        if (cur[d] != 0)
          {
            zero = false;
            break;
          }
      if (not zero)
        {
          result[idx] = cur;
          ++idx;
        }
      // Increment lex order (last axis fastest).
      for (std::size_t d = N; d-- > 0;)
        {
          if (cur[d] < static_cast<ca_index_t>(R))
            {
              ++cur[d];
              break;
            }
          cur[d] = -static_cast<ca_index_t>(R);
        }
    }
  return result;
}
 
template <std::size_t N, std::size_t R>
constexpr auto compute_von_neumann_offsets() noexcept
{
  constexpr std::size_t total = static_pow(2 * R + 1, N);
  constexpr std::size_t out_n = von_neumann_count<N, R>();
  std::array<Offset_Vec<N>, out_n> result{};
 
  Offset_Vec<N> cur{};
  for (std::size_t d = 0; d < N; ++d)
    cur[d] = -static_cast<ca_index_t>(R);
 
  std::size_t idx = 0;
  for (std::size_t i = 0; i < total; ++i)
    {
      std::size_t l1 = 0;
      for (std::size_t d = 0; d < N; ++d)
        l1 += static_cast<std::size_t>(cur[d] < 0 ? -cur[d] : cur[d]);
      if (l1 != 0 and l1 <= R)
        {
          result[idx] = cur;
          ++idx;
        }
      for (std::size_t d = N; d-- > 0;)
        {
          if (cur[d] < static_cast<ca_index_t>(R))
            {
              ++cur[d];
              break;
            }
          cur[d] = -static_cast<ca_index_t>(R);
        }
    }
  return result;
}
}  // namespace ca_neighborhood_detail
 
// -----------------------------------------------------------------------
// Moore neighborhood
// -----------------------------------------------------------------------
 
template <std::size_t N, std::size_t R = 1>
class Moore
{
  static_assert(N >= 1, "Moore requires N >= 1");
  static_assert(R >= 1, "Moore requires R >= 1");
 
public:
  static constexpr std::size_t rank_v = N;
  static constexpr std::size_t radius_v = R;
  static constexpr std::size_t size_v = ca_neighborhood_detail::static_pow(2 * R + 1, N) - 1;
  using coord_type = Coord_Vec<N>;
 
  static constexpr std::array<Offset_Vec<N>, size_v> offsets
    = ca_neighborhood_detail::compute_moore_offsets<N, R>();
 
  [[nodiscard]] static constexpr std::size_t radius() noexcept
  {
    return R;
  }
 
  [[nodiscard]] constexpr std::size_t size() const noexcept
  {
    return size_v;
  }
 
  template <typename F>
  constexpr void for_each_offset(const Coord_Vec<N> & center, F &&f) const
  {
    (void)center;
    for (std::size_t i = 0; i < size_v; ++i)
      f(offsets[i]);
  }
};
 
// -----------------------------------------------------------------------
// Von Neumann neighborhood
// -----------------------------------------------------------------------
 
template <std::size_t N, std::size_t R = 1>
class Von_Neumann
{
  static_assert(N >= 1, "Von_Neumann requires N >= 1");
  static_assert(R >= 1, "Von_Neumann requires R >= 1");
 
public:
  static constexpr std::size_t rank_v = N;
  static constexpr std::size_t radius_v = R;
  static constexpr std::size_t size_v = ca_neighborhood_detail::von_neumann_count<N, R>();
  using coord_type = Coord_Vec<N>;
 
  static constexpr std::array<Offset_Vec<N>, size_v> offsets
    = ca_neighborhood_detail::compute_von_neumann_offsets<N, R>();
 
  [[nodiscard]] static constexpr std::size_t radius() noexcept
  {
    return R;
  }
 
  [[nodiscard]] static constexpr std::size_t size() noexcept
  {
    return size_v;
  }
 
  template <typename F>
  constexpr void for_each_offset(const Coord_Vec<N> & center, F &&f) const
  {
    (void)center;
    for (std::size_t i = 0; i < size_v; ++i)
      f(offsets[i]);
  }
};
 
// -----------------------------------------------------------------------
// Hexagonal neighborhood (axial coordinates)
// -----------------------------------------------------------------------
 
class Hex_Neighborhood
{
public:
  static constexpr std::size_t rank_v = 2;
  static constexpr std::size_t radius_v = 1;
  static constexpr std::size_t size_v = 6;
  using coord_type = Coord_Vec<2>;
 
  static constexpr std::array<Offset_Vec<2>, 6> offsets
    = {{{1, 0}, {-1, 0}, {0, 1}, {0, -1}, {1, -1}, {-1, 1}}};
 
  [[nodiscard]] static constexpr std::size_t radius() noexcept
  {
    return 1;
  }
 
  [[nodiscard]] static constexpr std::size_t size() noexcept
  {
    return 6;
  }
 
  template <typename F>
  static constexpr void for_each_offset(const Coord_Vec<2> & center, F &&f)
  {
    (void)center;
    for (const auto &o : offsets)
      f(o);
  }
};
 
// -----------------------------------------------------------------------
// Triangular neighborhood
// -----------------------------------------------------------------------
 
class Triangular_Neighborhood
{
public:
  static constexpr std::size_t rank_v = 2;
  static constexpr std::size_t radius_v = 1;
  static constexpr std::size_t size_v = 3;
  using coord_type = Coord_Vec<2>;
 
  [[nodiscard]] static constexpr std::size_t radius() noexcept
  {
    return 1;
  }
 
  [[nodiscard]] static constexpr std::size_t size() noexcept
  {
    return 3;
  }
 
  template <typename F>
  static constexpr void for_each_offset(const Coord_Vec<2> &center, F &&f)
  {
    // Two parity-independent edge-neighbours along the j axis.
    f(Offset_Vec<2>{0, -1});
    f(Offset_Vec<2>{0, 1});
    // Third edge-neighbour points up or down depending on parity.
    const ca_index_t parity = (center[0] + center[1]) & 1;
    f(parity == 0 ? Offset_Vec<2>{1, 0} : Offset_Vec<2>{-1, 0});
  }
};
 
// -----------------------------------------------------------------------
// Custom neighborhood
// -----------------------------------------------------------------------
 
template <std::size_t N, std::size_t K, std::size_t Radius = 0>
class Custom_Neighborhood
{
  static_assert(N >= 1, "Custom_Neighborhood requires N >= 1");
  static_assert(K >= 1, "Custom_Neighborhood requires K >= 1");
 
  std::array<Offset_Vec<N>, K> offs_;
 
public:
  static constexpr std::size_t rank_v = N;
  static constexpr std::size_t size_v = K;
  static constexpr std::size_t radius_v = Radius;
  using coord_type = Coord_Vec<N>;
 
  constexpr explicit Custom_Neighborhood(std::array<Offset_Vec<N>, K> os) : offs_(os) {}
 
  constexpr explicit Custom_Neighborhood(const Offset_Vec<N> (&os)[K])
  {
    for (std::size_t i = 0; i < K; ++i)
      offs_[i] = os[i];
  }
 
  [[nodiscard]] constexpr std::size_t radius() const noexcept
  {
    std::size_t r = 0;
    for (std::size_t k = 0; k < K; ++k)
      for (std::size_t d = 0; d < N; ++d)
        {
          const auto v = offs_[k][d];
          const std::size_t a = static_cast<std::size_t>(v < 0 ? -v : v);
          if (a > r)
            r = a;
        }
    return r;
  }
 
  [[nodiscard]] constexpr std::size_t size() const noexcept
  {
    return K;
  }
 
  [[nodiscard]] constexpr const std::array<Offset_Vec<N>, K> &offsets() const noexcept
  {
    return offs_;
  }
 
  template <typename F>
  constexpr void for_each_offset(const Coord_Vec<N> & center, F &&f) const
  {
    (void)center;
    for (const auto &o : offs_)
      f(o);
  }
};
 
// -----------------------------------------------------------------------
// gather_neighbors: fills a caller-supplied buffer using at_safe.
// -----------------------------------------------------------------------
 
template <typename Nbh, typename L, typename T>
inline void gather_neighbors(const Nbh &nh,
                             const L &lat,
                             const typename L::coord_type &center,
                             std::span<T> out)
{
  ah_length_error_if(out.size() < nh.size()) << "gather_neighbors: output span has size "
                                             << out.size() << ", expected at least " << nh.size();
  using Coord = typename L::coord_type;
  std::size_t i = 0;
  nh.for_each_offset(center,
                     [&](const auto &off)
  {
    Coord c = center;
    for (std::size_t d = 0; d < L::rank; ++d)
      c[d] = static_cast<ca_index_t>(c[d] + off[d]);
    out[i] = static_cast<T>(lat.at_safe(c));
    ++i;
  });
}
 
}  // namespace CA
}  // namespace Aleph
 
#endif  // TPL_CA_NEIGHBORHOOD_H