tpl_ca_update_scheme.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_UPDATE_SCHEME_H
#define TPL_CA_UPDATE_SCHEME_H
 
#include <array>
#include <concepts>
#include <cstddef>
#include <cstdint>
#include <random>
#include <span>
#include <type_traits>
#include <utility>
 
#include <ah-errors.H>
#include <ca-rng.H>
#include <ca-traits.H>
#include <tpl_ca_concepts.H>
#include <tpl_ca_neighborhood.H>
 
namespace Aleph {
namespace CA {
 
namespace ca_scheme_detail {
 
template <typename O>
struct is_tile : std::false_type
{
};
 
template <std::size_t W, std::size_t H>
struct is_tile<Tile<W, H>> : std::true_type
{
};
 
template <typename O>
inline constexpr bool is_tile_v = is_tile<O>::value;
 
template <typename Coord, typename Cb>
inline void for_each_coord_1d(const ca_size_t n0, Cb &&cb)
{
  Coord c{};
  for (ca_index_t i = 0; i < static_cast<ca_index_t>(n0); ++i)
    {
      c[0] = i;
      cb(c);
    }
}
 
template <typename Coord, typename Cb>
inline void for_each_coord_2d_row(const ca_size_t n0, const ca_size_t n1, Cb &&cb)
{
  Coord c{};
  for (ca_index_t i = 0; i < static_cast<ca_index_t>(n0); ++i)
    for (ca_index_t j = 0; j < static_cast<ca_index_t>(n1); ++j)
      {
        c[0] = i;
        c[1] = j;
        cb(c);
      }
}
 
template <typename Coord, typename Cb>
inline void for_each_coord_2d_col(const ca_size_t n0, const ca_size_t n1, Cb &&cb)
{
  Coord c{};
  for (ca_index_t j = 0; j < static_cast<ca_index_t>(n1); ++j)
    for (ca_index_t i = 0; i < static_cast<ca_index_t>(n0); ++i)
      {
        c[0] = i;
        c[1] = j;
        cb(c);
      }
}
 
template <std::size_t W, std::size_t H, typename Coord, typename Cb>
inline void for_each_coord_2d_tile(const ca_size_t n0, const ca_size_t n1, Cb &&cb)
{
  Coord c{};
  for (ca_size_t ti = 0; ti < n0; ti += H)
    for (ca_size_t tj = 0; tj < n1; tj += W)
      {
        const ca_size_t i_end = ti + H < n0 ? ti + H : n0;
        const ca_size_t j_end = tj + W < n1 ? tj + W : n1;
        for (ca_size_t i = ti; i < i_end; ++i)
          for (ca_size_t j = tj; j < j_end; ++j)
            {
              c[0] = static_cast<ca_index_t>(i);
              c[1] = static_cast<ca_index_t>(j);
              cb(c);
            }
      }
}
 
template <typename Coord, typename Cb>
inline void for_each_coord_3d_row(const ca_size_t n0, const ca_size_t n1, const ca_size_t n2, Cb &&cb)
{
  Coord c{};
  for (ca_index_t i = 0; i < static_cast<ca_index_t>(n0); ++i)
    for (ca_index_t j = 0; j < static_cast<ca_index_t>(n1); ++j)
      for (ca_index_t k = 0; k < static_cast<ca_index_t>(n2); ++k)
        {
          c[0] = i;
          c[1] = j;
          c[2] = k;
          cb(c);
        }
}
 
template <typename Order, typename Lattice, typename Cb>
inline void for_each_coord(const Lattice &lat, Cb &&cb)
{
  using Coord = typename Lattice::coord_type;
  if constexpr (Lattice::rank == 1)
    {
      for_each_coord_1d<Coord>(lat.size(0), std::forward<Cb>(cb));
    }
  else if constexpr (Lattice::rank == 2)
    {
      if constexpr (is_tile_v<Order>)
        for_each_coord_2d_tile<Order::width, Order::height, Coord>(lat.size(0), lat.size(1),
                                                                   std::forward<Cb>(cb));
      else if constexpr (std::is_same_v<Order, ColumnMajor>)
        for_each_coord_2d_col<Coord>(lat.size(0), lat.size(1), std::forward<Cb>(cb));
      else
        for_each_coord_2d_row<Coord>(lat.size(0), lat.size(1), std::forward<Cb>(cb));
    }
  else if constexpr (Lattice::rank == 3)
    {
      for_each_coord_3d_row<Coord>(lat.size(0), lat.size(1), lat.size(2), std::forward<Cb>(cb));
    }
  else
    {
      static_assert(Lattice::rank <= 3, "Update schemes currently support rank <= 3");
    }
}
 
template <typename Lattice>
inline ca_size_t cell_count(const Lattice &lat) noexcept
{
  ca_size_t n = 1;
  for (std::size_t d = 0; d < Lattice::rank; ++d)
    n *= lat.size(d);
  return n;
}
 
template <std::size_t Rank>
inline Coord_Vec<Rank> coord_from_hash(std::uint64_t h, const std::array<ca_size_t, Rank> &ext) noexcept
{
  Coord_Vec<Rank> c{};
  for (std::size_t d = 0; d < Rank; ++d)
    {
      const ca_size_t n = ext[d] == 0 ? 1 : ext[d];
      c[d] = static_cast<ca_index_t>(h % n);
      h /= n == 0 ? 1 : n;
      h = splitmix64(h ^ (static_cast<std::uint64_t>(d) + 1));
    }
  return c;
}
 
}  // namespace ca_scheme_detail
 
// =======================================================================
// Synchronous_Update
// =======================================================================
 
template <typename Order = RowMajor>
struct Synchronous_Update
{
  static constexpr bool requires_double_buffer = true;
  static constexpr bool requires_block_rule = false;
 
  template <typename Engine, typename Rule>
  void apply(Engine &e, Rule &r) const
  {
    using Lattice = typename Engine::lattice_type;
    using State = typename Lattice::state_type;
    using Coord = typename Lattice::coord_type;
    using Neighborhood = typename Engine::neighborhood_type;
    constexpr std::size_t K = Neighborhood::size_v;
 
    auto &cur = e.current_buffer();
    auto &nxt = e.next_buffer();
    const auto &nh = e.neighborhood();
    const std::size_t step = e.step_count();
 
    std::array<State, K == 0 ? 1 : K> nbuf{};
    const std::size_t neigh_sz = (std::min) (nh.size(), static_cast<std::size_t>(nbuf.size()));
 
    ca_scheme_detail::for_each_coord<Order>(cur, [&](const Coord &c)
    {
      gather_neighbors(nh, cur, c, std::span<State>(nbuf.data(), neigh_sz));
      const Cell_Context<Lattice::rank> ctx{step, c};
      nxt.set(c, apply_rule(r, cur.at(c), Neighbor_View<State>(nbuf.data(), neigh_sz), ctx));
    });
 
    e.swap_buffers();
  }
};
 
// =======================================================================
// Sequential_Update
// =======================================================================
 
template <typename Order = RowMajor>
struct Sequential_Update
{
  static constexpr bool requires_double_buffer = false;
  static constexpr bool requires_block_rule = false;
 
  template <typename Engine, typename Rule>
  void apply(Engine &e, Rule &r) const
  {
    using Lattice = typename Engine::lattice_type;
    using State = typename Lattice::state_type;
    using Coord = typename Lattice::coord_type;
    using Neighborhood = typename Engine::neighborhood_type;
    constexpr std::size_t K = Neighborhood::size_v;
 
    auto &cur = e.current_buffer();
    const auto &nh = e.neighborhood();
    const std::size_t step = e.step_count();
 
    std::array<State, K == 0 ? 1 : K> nbuf{};
    const std::size_t neigh_sz = (std::min) (nh.size(), static_cast<std::size_t>(nbuf.size()));
 
    ca_scheme_detail::for_each_coord<Order>(cur, [&](const Coord &c)
    {
      gather_neighbors(nh, cur, c, std::span<State>(nbuf.data(), neigh_sz));
      const Cell_Context<Lattice::rank> ctx{step, c};
      const State next = apply_rule(r, cur.at(c), Neighbor_View<State>(nbuf.data(), neigh_sz), ctx);
      cur.set(c, next);
    });
  }
};
 
// =======================================================================
// Random_Asynchronous_Update
// =======================================================================
 
template <typename Engine_Type = std::mt19937_64>
struct Random_Asynchronous_Update
{
  static constexpr bool requires_double_buffer = false;
  static constexpr bool requires_block_rule = false;
 
  std::uint64_t master_seed = 0;
  std::size_t sub_steps_per_step = 0;
 
  constexpr Random_Asynchronous_Update(const std::uint64_t seed = 0,
                                       const std::size_t sub_steps = 0) noexcept
    : master_seed(seed), sub_steps_per_step(sub_steps)
  {}
 
  template <typename Engine, typename Rule>
  void apply(Engine &e, Rule &r) const
  {
    using Lattice = typename Engine::lattice_type;
    using State = typename Lattice::state_type;
    using Coord = typename Lattice::coord_type;
    using Neighborhood = typename Engine::neighborhood_type;
    constexpr std::size_t K = Neighborhood::size_v;
 
    auto &cur = e.current_buffer();
    const auto &nh = e.neighborhood();
    const std::size_t step = e.step_count();
 
    const ca_size_t n_cells = ca_scheme_detail::cell_count(cur);
    const std::size_t subs
      = sub_steps_per_step == 0 ? static_cast<std::size_t>(n_cells) : sub_steps_per_step;
 
    std::array<State, K == 0 ? 1 : K> nbuf{};
    const std::size_t neigh_sz = (std::min) (nh.size(), static_cast<std::size_t>(nbuf.size()));
 
    for (std::size_t sub = 0; sub < subs; ++sub)
      {
        // Stable per-(step, sub) selection seed independent of the
        // rule's own per-cell sub-stream.
        const std::uint64_t pick_seed = mix_seed(
          mix_seed(master_seed, static_cast<std::uint64_t>(step)), static_cast<std::uint64_t>(sub));
        Coord c = ca_scheme_detail::coord_from_hash<Lattice::rank>(pick_seed, cur.extents());
 
        gather_neighbors(nh, cur, c, std::span<State>(nbuf.data(), neigh_sz));
        const Cell_Context<Lattice::rank> ctx{step, c};
        const State next = apply_rule(r, cur.at(c), Neighbor_View<State>(nbuf.data(), neigh_sz), ctx);
        cur.set(c, next);
      }
  }
};
 
// =======================================================================
// Block_Synchronous_Update
// =======================================================================
 
template <std::size_t BlockSize = 4, typename Order = RowMajor>
struct Block_Synchronous_Update
{
  static_assert(BlockSize >= 1, "Block_Synchronous_Update requires BlockSize >= 1");
 
  static constexpr bool requires_double_buffer = true;
  static constexpr bool requires_block_rule = false;
 
  static constexpr std::size_t block_size = BlockSize;
 
  template <typename Engine, typename Rule>
  void apply(Engine &e, Rule &r) const
  {
    using Lattice = typename Engine::lattice_type;
    using State = typename Lattice::state_type;
    using Coord = typename Lattice::coord_type;
    using Neighborhood = typename Engine::neighborhood_type;
    constexpr std::size_t K = Neighborhood::size_v;
 
    auto &cur = e.current_buffer();
    auto &nxt = e.next_buffer();
    const auto &nh = e.neighborhood();
    const std::size_t step = e.step_count();
 
    // Number of blocks along each axis.
    std::array<ca_size_t, Lattice::rank> blocks_per_axis{};
    ca_size_t total_blocks = 1;
    for (std::size_t d = 0; d < Lattice::rank; ++d)
      {
        const ca_size_t n = cur.size(d);
        blocks_per_axis[d] = (n + BlockSize - 1) / BlockSize;
        total_blocks *= (blocks_per_axis[d] == 0 ? 1 : blocks_per_axis[d]);
      }
 
    if (total_blocks == 0)
      total_blocks = 1;
    const std::size_t active = step % total_blocks;
 
    // Decompose the active block index into per-axis tile coordinates.
    std::array<ca_size_t, Lattice::rank> tile_coord{};
    {
      std::size_t rem = active;
      for (std::size_t d = Lattice::rank; d-- > 0;)
        {
          const ca_size_t bn = blocks_per_axis[d] == 0 ? 1 : blocks_per_axis[d];
          tile_coord[d] = rem % bn;
          rem /= bn;
        }
    }
 
    // First, copy current → next so cells outside the active block
    // carry over unchanged after the swap.
    ca_scheme_detail::for_each_coord<Order>(cur, [&](const Coord &c)
    {
      nxt.set(c, cur.at(c));
    });
 
    std::array<State, K == 0 ? 1 : K> nbuf{};
    const std::size_t neigh_sz = (std::min) (nh.size(), static_cast<std::size_t>(nbuf.size()));
 
    // Iterate over the active block only.
    if constexpr (Lattice::rank == 1)
      {
        const ca_size_t i0 = tile_coord[0] * BlockSize;
        const ca_size_t i_end = (std::min) (i0 + BlockSize, cur.size(0));
        for (ca_size_t i = i0; i < i_end; ++i)
          {
            Coord c{static_cast<ca_index_t>(i)};
            gather_neighbors(nh, cur, c, std::span<State>(nbuf.data(), neigh_sz));
            const Cell_Context<Lattice::rank> ctx{step, c};
            nxt.set(c, apply_rule(r, cur.at(c), Neighbor_View<State>(nbuf.data(), neigh_sz), ctx));
          }
      }
    else if constexpr (Lattice::rank == 2)
      {
        const ca_size_t i0 = tile_coord[0] * BlockSize;
        const ca_size_t j0 = tile_coord[1] * BlockSize;
        const ca_size_t i_end = (std::min) (i0 + BlockSize, cur.size(0));
        const ca_size_t j_end = (std::min) (j0 + BlockSize, cur.size(1));
        for (ca_size_t i = i0; i < i_end; ++i)
          for (ca_size_t j = j0; j < j_end; ++j)
            {
              Coord c{static_cast<ca_index_t>(i), static_cast<ca_index_t>(j)};
              gather_neighbors(nh, cur, c, std::span<State>(nbuf.data(), neigh_sz));
              const Cell_Context<Lattice::rank> ctx{step, c};
              nxt.set(c, apply_rule(r, cur.at(c), Neighbor_View<State>(nbuf.data(), neigh_sz), ctx));
            }
      }
    else if constexpr (Lattice::rank == 3)
      {
        const ca_size_t i0 = tile_coord[0] * BlockSize;
        const ca_size_t j0 = tile_coord[1] * BlockSize;
        const ca_size_t k0 = tile_coord[2] * BlockSize;
        const ca_size_t i_end = (std::min) (i0 + BlockSize, cur.size(0));
        const ca_size_t j_end = (std::min) (j0 + BlockSize, cur.size(1));
        const ca_size_t k_end = (std::min) (k0 + BlockSize, cur.size(2));
        for (ca_size_t i = i0; i < i_end; ++i)
          for (ca_size_t j = j0; j < j_end; ++j)
            for (ca_size_t k = k0; k < k_end; ++k)
              {
                Coord c{static_cast<ca_index_t>(i), static_cast<ca_index_t>(j),
                        static_cast<ca_index_t>(k)};
                gather_neighbors(nh, cur, c, std::span<State>(nbuf.data(), neigh_sz));
                const Cell_Context<Lattice::rank> ctx{step, c};
                nxt.set(c,
                        apply_rule(r, cur.at(c), Neighbor_View<State>(nbuf.data(), neigh_sz), ctx));
              }
      }
 
    e.swap_buffers();
  }
};
 
// =======================================================================
// Margolus_Update
// =======================================================================
 
struct Margolus_Update
{
  static constexpr bool requires_double_buffer = false;
  static constexpr bool requires_block_rule = true;
 
  [[nodiscard]] static constexpr ca_index_t origin_for_step(const std::size_t k) noexcept
  {
    return (k % 2 == 0) ? 0 : 1;
  }
 
  template <typename Engine, typename Block_Rule>
  void apply_with_origin(Engine &e, Block_Rule &r, const ca_index_t origin) const
  {
    using Lattice = typename Engine::lattice_type;
    using State = typename Lattice::state_type;
    using Coord = typename Lattice::coord_type;
    static_assert(Lattice::rank == 2, "Margolus_Update requires a 2D lattice");
 
    auto &lat = e.current_buffer();
    const ca_size_t n0 = lat.size(0);
    const ca_size_t n1 = lat.size(1);
    if (n0 < 2 or n1 < 2)
      return;
 
    // Sweep every 2×2 block of the active partition.
    for (ca_index_t i = origin; i + 1 < static_cast<ca_index_t>(n0); i += 2)
      for (ca_index_t j = origin; j + 1 < static_cast<ca_index_t>(n1); j += 2)
        {
          const Coord nw{i, j};
          const Coord ne{i, j + 1};
          const Coord sw{i + 1, j};
          const Coord se{i + 1, j + 1};
 
          std::array<State, 4> in{lat.at(nw), lat.at(ne), lat.at(sw), lat.at(se)};
          const std::array<State, 4> out = r(in);
 
          lat.set(nw, out[0]);
          lat.set(ne, out[1]);
          lat.set(sw, out[2]);
          lat.set(se, out[3]);
        }
  }
 
  template <typename Engine, typename Block_Rule>
  void apply(Engine &e, Block_Rule &r) const
  {
    apply_with_origin(e, r, origin_for_step(e.step_count()));
  }
};
 
}  // namespace CA
}  // namespace Aleph
 
#endif  // TPL_CA_UPDATE_SCHEME_H