diff --git a/demos/FiniteVolume/linear_convection_obstacle.cpp b/demos/FiniteVolume/linear_convection_obstacle.cpp index a39c513e0..f8329ed83 100644 --- a/demos/FiniteVolume/linear_convection_obstacle.cpp +++ b/demos/FiniteVolume/linear_convection_obstacle.cpp @@ -136,7 +136,8 @@ int main(int argc, char* argv[]) } double t = 0; - while (t != Tf) + //~ while (t != Tf) + for (int ite = 0; ite != 23; ++ite) { // Move to next timestep t += dt; diff --git a/demos/tutorial/set_operator.cpp b/demos/tutorial/set_operator.cpp index 4f46700b1..82a4bd53e 100644 --- a/demos/tutorial/set_operator.cpp +++ b/demos/tutorial/set_operator.cpp @@ -98,7 +98,7 @@ int main() u[cell] = cell.indices[0]; }); - auto subset1 = samurai::intersection(ca[0], samurai::contraction(ca[1], 1)); + auto subset1 = samurai::intersection(ca[0], samurai::contract(ca[1], 1)); subset1.on(0)( [&](const auto& i, auto) { diff --git a/include/samurai/algorithm.hpp b/include/samurai/algorithm.hpp index bd6c21617..9c262b4a4 100644 --- a/include/samurai/algorithm.hpp +++ b/include/samurai/algorithm.hpp @@ -14,6 +14,9 @@ #include "cell.hpp" #include "mesh_holder.hpp" +#include "concepts.hpp" +#include "subset/node.hpp" + using namespace xt::placeholders; namespace samurai @@ -111,18 +114,28 @@ namespace samurai } } - template + template inline void for_each_interval(const Mesh& mesh, Func&& f) { using mesh_id_t = typename Mesh::config::mesh_id_t; for_each_interval(mesh[mesh_id_t::cells], std::forward(f)); } - template - class Subset; - - template - inline void for_each_interval(Subset& set, Func&& f) + //~ template + //~ class Subset; + + //~ template + //~ inline void for_each_interval(Subset& set, Func&& f) + //~ { + //~ set( + //~ [&](const auto& i, const auto& index) + //~ { + //~ f(set.level(), i, index); + //~ }); + //~ } + + template + inline void for_each_interval(const SetBase& set, Func&& f) { set( [&](const auto& i, const auto& index) diff --git a/include/samurai/boundary.hpp b/include/samurai/boundary.hpp index ac505f65e..b22d62ae2 100644 --- a/include/samurai/boundary.hpp +++ b/include/samurai/boundary.hpp @@ -9,7 +9,7 @@ namespace samurai { using mesh_id_t = typename Mesh::mesh_id_t; - auto& cells = mesh[mesh_id_t::cells][level]; + const auto& cells = mesh[mesh_id_t::cells][level]; return difference(cells, translate(self(domain).on(level), -layer_width * direction)); } @@ -43,7 +43,7 @@ namespace samurai { using mesh_id_t = typename Mesh::mesh_id_t; - auto& cells = mesh[mesh_id_t::cells][level]; + const auto& cells = mesh[mesh_id_t::cells][level]; return difference(cells, contract(self(mesh.domain()).on(level), 1)); } diff --git a/include/samurai/interval.hpp b/include/samurai/interval.hpp index 97ef0cd5d..c84a60263 100644 --- a/include/samurai/interval.hpp +++ b/include/samurai/interval.hpp @@ -362,7 +362,8 @@ namespace samurai template inline bool operator==(const Interval& i1, const Interval& i2) { - return !(i1.start != i2.start || i1.end != i2.end || i1.step != i2.step || i1.index != i2.index); + //~ return !(i1.start != i2.start || i1.end != i2.end || i1.step != i2.step || i1.index != i2.index); + return !(i1.start != i2.start || i1.end != i2.end || i1.step != i2.step); } template diff --git a/include/samurai/level_cell_array.hpp b/include/samurai/level_cell_array.hpp index 48b79c717..18b8c1eba 100644 --- a/include/samurai/level_cell_array.hpp +++ b/include/samurai/level_cell_array.hpp @@ -94,8 +94,11 @@ namespace samurai LevelCellArray() = default; LevelCellArray(const LevelCellList& lcl); - template - LevelCellArray(Subset set); + template + explicit LevelCellArray(const SetBase& set); + + template + LevelCellArray(const SetBase& set, const coords_t& origin_point, const double scaling_factor); LevelCellArray(std::size_t level, const Box& box); LevelCellArray(std::size_t level, @@ -344,9 +347,21 @@ namespace samurai } } + //~ template + //~ template + //~ inline LevelCellArray::LevelCellArray(Subset set) + //~ : m_level(set.level()) + //~ { + //~ set( + //~ [this](const auto& i, const auto& index) + //~ { + //~ add_interval_back(i, index); + //~ }); + //~ } + template - template - inline LevelCellArray::LevelCellArray(Subset set) + template + inline LevelCellArray::LevelCellArray(const SetBase& set) : m_level(set.level()) { set( @@ -356,6 +371,32 @@ namespace samurai }); } + template + template + inline LevelCellArray::LevelCellArray(const SetBase& set, const coords_t& origin_point, const double scaling_factor) + : m_level(set.level()) + , m_origin_point(origin_point) + , m_scaling_factor(scaling_factor) + { + set( + [this](const auto& i, const auto& index) + { + add_interval_back(i, index); + }); + } + + //~ template + //~ template + //~ inline LevelCellArray::LevelCellArray(Subset set) + //~ : m_level(set.level()) + //~ { + //~ set( + //~ [this](const auto& i, const auto& index) + //~ { + //~ add_interval_back(i, index); + //~ }); + //~ } + template inline LevelCellArray::LevelCellArray(std::size_t level, const Box& box) : m_level{level} diff --git a/include/samurai/list_of_intervals.hpp b/include/samurai/list_of_intervals.hpp index d3889f847..ee8d90e12 100644 --- a/include/samurai/list_of_intervals.hpp +++ b/include/samurai/list_of_intervals.hpp @@ -62,6 +62,11 @@ namespace samurai void add_point(value_t point); void add_interval(const interval_t& interval); + + void clear() + { + std::forward_list>::clear(); + } }; //////////////////////////////////// diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index 1fa2192f7..ed17be689 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -815,16 +815,18 @@ namespace samurai { if constexpr (dim == 2) { - m_corners.push_back(difference( - m_domain, - union_(translate(m_domain, direction_t{-direction[0], 0}), translate(m_domain, direction_t{0, -direction[1]})))); + m_corners.push_back(difference(m_domain, + union_(translate(m_domain, direction_t{-direction[0], 0}), + translate(m_domain, direction_t{0, -direction[1]}))) + .to_lca()); } else if constexpr (dim == 3) { m_corners.push_back(difference(m_domain, union_(translate(m_domain, direction_t{-direction[0], 0, 0}), translate(m_domain, direction_t{0, -direction[1], 0}), - translate(m_domain, direction_t{0, 0, -direction[2]})))); + translate(m_domain, direction_t{0, 0, -direction[2]}))) + .to_lca()); } }); } diff --git a/include/samurai/samurai_config.hpp b/include/samurai/samurai_config.hpp index b4394f61a..f5ac59a5a 100644 --- a/include/samurai/samurai_config.hpp +++ b/include/samurai/samurai_config.hpp @@ -20,6 +20,8 @@ namespace samurai static constexpr std::size_t graduation_width = 1; static constexpr std::size_t prediction_order = 1; + static constexpr bool prediction_with_list_of_intervals = false; + using index_t = signed long long int; using value_t = int; using interval_t = Interval; diff --git a/include/samurai/subset/apply.hpp b/include/samurai/subset/apply.hpp index b8652afe2..c75ebfb62 100644 --- a/include/samurai/subset/apply.hpp +++ b/include/samurai/subset/apply.hpp @@ -3,119 +3,50 @@ #pragma once -#include "concepts.hpp" -#include "utils.hpp" +#include "set_base.hpp" namespace samurai { namespace detail { - template - bool apply_impl(Set&& global_set, Func&& func, Container& index) + template + void apply_rec(const SetBase& set, Func&& func, Index& index, std::integral_constant d_ic) { - auto set = global_set.template get_local_set(global_set.level(), index); - auto start_and_stop = global_set.template get_start_and_stop_function(); + using traverser_t = typename Set::template traverser_t; + using current_interval_t = typename traverser_t::current_interval_t; - if constexpr (dim != 1) + set.init_get_traverser_work(1, d_ic); + + for (traverser_t traverser = set.get_traverser(index, d_ic); !traverser.is_empty(); traverser.next_interval()) { - auto func_int = [&](const auto& interval) + current_interval_t interval = traverser.current_interval(); + + if constexpr (d == 0) { - for (auto i = interval.start; i < interval.end; ++i) + func(interval, index); + } + else + { + for (index[d - 1] = interval.start; index[d - 1] != interval.end; ++index[d - 1]) { - index[dim - 2] = i; - if (apply_impl(std::forward(global_set), std::forward(func), index)) - { - return true; - } + apply_rec(set, std::forward(func), index, std::integral_constant{}); } - return false; - }; - return apply(set, start_and_stop, func_int); - } - else - { - auto func_int = [&](const auto& interval) - { - return func(interval, index); - }; - return apply(set, start_and_stop, func_int); + } } + + set.clear_get_traverser_work(d_ic); } } template - void apply(Set&& global_set, Func&& user_func) + void apply(const SetBase& set, Func&& func) { - constexpr std::size_t dim = std::decay_t::dim; - xt::xtensor_fixed> index; - - auto func = [&](const auto& interval, const auto& yz) - { - user_func(interval, yz); - return false; - }; + constexpr std::size_t dim = Set::dim; - if (global_set.exist()) - { - detail::apply_impl(std::forward(global_set), func, index); - } - } - - template - bool empty_check(Set&& global_set) - { - constexpr std::size_t dim = std::decay_t::dim; xt::xtensor_fixed> index; - - auto func = [](const auto&, const auto&) + if (set.exist()) { - return true; - }; - - if (global_set.exist()) - { - return !detail::apply_impl(std::forward(global_set), func, index); - } - return true; - } - - template - requires IsSetOp || IsIntervalListVisitor - bool apply(Set&& set, StartEnd&& start_and_stop, Func&& func) - { - using interval_t = typename std::decay_t::interval_t; - using value_t = typename interval_t::value_t; - - interval_t result; - int r_ipos = 0; - set.next(0, std::forward(start_and_stop)); - auto scan = set.min(); - - while (scan < sentinel && !set.is_empty()) - { - bool is_in = set.is_in(scan); - - if (is_in && r_ipos == 0) - { - result.start = scan; - r_ipos = 1; - } - else if (!is_in && r_ipos == 1) - { - result.end = scan; - r_ipos = 0; - - auto true_result = set.shift() >= 0 ? result >> static_cast(set.shift()) - : result << -static_cast(set.shift()); - if (func(true_result)) - { - return true; - } - } - - set.next(scan, std::forward(start_and_stop)); - scan = set.min(); + detail::apply_rec(set, std::forward(func), index, std::integral_constant{}); } - return false; } } diff --git a/include/samurai/subset/box_view.hpp b/include/samurai/subset/box_view.hpp new file mode 100644 index 000000000..b2b6153b3 --- /dev/null +++ b/include/samurai/subset/box_view.hpp @@ -0,0 +1,83 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_base.hpp" +#include "traversers/box_traverser.hpp" + +namespace samurai +{ + + template + class BoxView; + + template + struct SetTraits> + { + static_assert(std::same_as, B>); + + template + using traverser_t = BoxTraverser; + + static constexpr std::size_t dim() + { + return B::dim; + } + }; + + template + class BoxView : public SetBase> + { + using Self = BoxView; + + public: + + SAMURAI_SET_TYPEDEFS + + BoxView(const std::size_t level, const B& box) + : m_level(level) + , m_box(box) + { + } + + inline std::size_t level_impl() const + { + return m_level; + } + + inline bool exist_impl() const + { + return m_box.is_valid(); + } + + inline bool empty_impl() const + { + return !exist_impl(); + } + + template + inline traverser_t get_traverser_impl(const index_t& index, std::integral_constant) const + { + return (m_box.min_corner()[d + 1] <= index[d] && index[d] < m_box.max_corner()[d + 1]) + ? traverser_t(m_box.min_corner()[d], m_box.max_corner()[d]) + : traverser_t(0, 0); + } + + template + constexpr inline void init_get_traverser_work_impl(const std::size_t, std::integral_constant) const + { + } + + template + inline void clear_get_traverser_work_impl(std::integral_constant) const + { + } + + private: + + std::size_t m_level; + const B& m_box; + }; + +} // namespace samurai diff --git a/include/samurai/subset/concepts.hpp b/include/samurai/subset/concepts.hpp deleted file mode 100644 index 4fcdc6125..000000000 --- a/include/samurai/subset/concepts.hpp +++ /dev/null @@ -1,44 +0,0 @@ -// Copyright 2018-2025 the samurai's authors -// SPDX-License-Identifier: BSD-3-Clause - -#pragma once - -#include -#include - -namespace samurai -{ - template - class LevelCellArray; - // } - - // namespace samurai::experimental - // { - template - class SetTraverser; - - template - class IntervalListVisitor; - - template - struct is_setop : std::false_type - { - }; - - template - struct is_setop> : std::true_type - { - }; - - template - constexpr bool is_setop_v{is_setop>::value}; - - template - concept IsSetOp = is_setop_v; - - template - concept IsIntervalListVisitor = std::is_base_of_v::container_t>, std::decay_t>; - - template - concept IsLCA = std::same_as, T>; -} diff --git a/include/samurai/subset/contraction.hpp b/include/samurai/subset/contraction.hpp new file mode 100644 index 000000000..784e66144 --- /dev/null +++ b/include/samurai/subset/contraction.hpp @@ -0,0 +1,130 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_base.hpp" +#include "traversers/contraction_traverser.hpp" + +namespace samurai +{ + + template + class Contraction; + + template + struct SetTraits> + { + static_assert(IsSet::value); + + template + using traverser_t = ContractionTraverser>; + + static constexpr std::size_t dim() + { + return Set::dim; + } + }; + + template + class Contraction : public SetBase> + { + using Self = Contraction; + + public: + + SAMURAI_SET_TYPEDEFS + + using contraction_t = std::array; + using do_contraction_t = std::array; + + Contraction(const Set& set, const contraction_t& contraction) + : m_set(set) + , m_contraction(contraction) + { + assert(std::all_of(m_contraction.cbegin(), + m_contraction.cend(), + [](const contraction_t& c) + { + return c >= 0; + })); + } + + Contraction(const Set& set, const value_t contraction) + : m_set(set) + { + assert(contraction >= 0); + std::fill(m_contraction.begin(), m_contraction.end(), contraction); + } + + Contraction(const Set& set, const value_t contraction, const do_contraction_t& do_contraction) + : m_set(set) + { + for (std::size_t i = 0; i != m_contraction.size(); ++i) + { + m_contraction[i] = contraction * do_contraction[i]; + } + } + + inline std::size_t level_impl() const + { + return m_set.level(); + } + + inline bool exist_impl() const + { + return m_set.exist(); + } + + inline bool empty_impl() const + { + return Base::empty_default_impl(); + } + + template + inline void init_get_traverser_work_impl(const std::size_t n_traversers, std::integral_constant d_ic) const + { + m_set.init_get_traverser_work(n_traversers, d_ic); + } + + template + inline void clear_get_traverser_work_impl(std::integral_constant d_ic) const + { + m_set.clear_get_traverser_work(d_ic); + } + + template + inline traverser_t get_traverser_impl(const index_t& index, std::integral_constant d_ic) const + { + return traverser_t(m_set.get_traverser(index, d_ic), m_contraction[d]); + } + + private: + + Set m_set; + contraction_t m_contraction; + }; + + template + auto contract(const Set& set, const typename Contraction>::contraction_t& contraction) + { + return Contraction(self(set), contraction); + } + + template + auto contract(const Set& set, const typename Contraction>::value_t& contraction) + { + return Contraction(self(set), contraction); + } + + template + auto contract(const Set& set, + const typename Contraction>::value_t& contraction, + const typename Contraction>::do_contraction_t& do_contraction) // idk how to make this + // more readable, + // perhaps a traits... + { + return Contraction(self(set), contraction, do_contraction); + } + +} // namespace samurai diff --git a/include/samurai/subset/expansion.hpp b/include/samurai/subset/expansion.hpp new file mode 100644 index 000000000..d873036b4 --- /dev/null +++ b/include/samurai/subset/expansion.hpp @@ -0,0 +1,211 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_base.hpp" +#include "traversers/expansion_traverser.hpp" +#include "traversers/last_dim_expansion_traverser.hpp" + +namespace samurai +{ + + template + class Expansion; + + template + struct SetTraits> + { + static_assert(IsSet::value); + + template + using traverser_t = std::conditional_t>, + ExpansionTraverser>>; + + static constexpr std::size_t dim() + { + return Set::dim; + } + }; + + namespace detail + { + template + struct ExpansionWork; + + template + struct ExpansionWork> + { + template + using child_traverser_t = typename Set::template traverser_t; + + template + using array_of_child_traverser_t = std::vector>; + + using Type = std::tuple...>; + + static_assert(std::tuple_size::value == Set::dim - 1); + }; + + } // namespace detail + + template + class Expansion : public SetBase> + { + using Self = Expansion; + using ChildTraverserArray = detail::ExpansionWork>::Type; + + public: + + SAMURAI_SET_TYPEDEFS + + using expansion_t = std::array; + using do_expansion_t = std::array; + + Expansion(const Set& set, const expansion_t& expansions) + : m_set(set) + , m_expansions(expansions) + { + } + + Expansion(const Set& set, const value_t expansion) + : m_set(set) + { + m_expansions.fill(expansion); + } + + Expansion(const Set& set, const value_t expansion, const do_expansion_t& do_expansion) + : m_set(set) + { + for (std::size_t i = 0; i != m_expansions.size(); ++i) + { + m_expansions[i] = expansion * do_expansion[i]; + } + } + + // we need to define a custom copy and move constructor because + // we do not want to copy m_work_offsetRanges + Expansion(const Expansion& other) + : m_set(other.m_set) + , m_expansions(other.m_expansions) + { + } + + Expansion(Expansion&& other) + : m_set(other.m_set) + , m_expansions(other.m_expansions) + { + } + + inline std::size_t level_impl() const + { + return m_set.level(); + } + + inline bool exist_impl() const + { + return m_set.exist(); + } + + inline bool empty_impl() const + { + return m_set.empty(); + } + + template + inline void init_get_traverser_work_impl(const std::size_t n_traversers, std::integral_constant d_ic) const + { + if constexpr (d == Base::dim - 1) + { + m_set.init_get_traverser_work_impl(n_traversers, d_ic); + } + else + { + const std::size_t my_work_size = n_traversers * 2 * std::size_t(m_expansions[d + 1] + 1); + + auto& childTraversers = std::get(m_work_childTraversers); + childTraversers.clear(); + childTraversers.reserve(my_work_size); + + m_set.init_get_traverser_work_impl(my_work_size, d_ic); + } + } + + template + inline void clear_get_traverser_work_impl(std::integral_constant d_ic) const + { + if constexpr (d != Base::dim - 1) + { + auto& childTraversers = std::get(m_work_childTraversers); + + childTraversers.clear(); + } + m_set.clear_get_traverser_work(d_ic); + } + + template + inline traverser_t get_traverser_impl(const index_t& index, std::integral_constant d_ic) const + { + if constexpr (d == Base::dim - 1) + { + return traverser_t(m_set.get_traverser(index, d_ic), m_expansions[d]); + } + else + { + auto& childTraversers = std::get(m_work_childTraversers); + + xt::xtensor_fixed> tmp_index(index); + + const auto childTraversers_begin = childTraversers.end(); + + for (value_t width = 0; width != m_expansions[d + 1] + 1; ++width) + { + tmp_index[d + 1] = index[d + 1] + width; + childTraversers.push_back(m_set.get_traverser(tmp_index, d_ic)); + if (childTraversers.back().is_empty()) + { + childTraversers.pop_back(); + } + + tmp_index[d + 1] = index[d + 1] - width; + childTraversers.push_back(m_set.get_traverser(tmp_index, d_ic)); + if (childTraversers.back().is_empty()) + { + childTraversers.pop_back(); + } + } + + return traverser_t(childTraversers_begin, childTraversers.end(), m_expansions[d]); + } + } + + private: + + Set m_set; + expansion_t m_expansions; + + mutable ChildTraverserArray m_work_childTraversers; + }; + + template + auto expand(const Set& set, const typename Contraction>::contraction_t& expansions) + { + return Expansion(self(set), expansions); + } + + template + auto expand(const Set& set, const typename Contraction>::value_t& expansion) + { + return Expansion(self(set), expansion); + } + + template + auto expand(const Set& set, + const typename Contraction>::value_t& expansion, + const typename Contraction>::do_expansion_t& do_expansion) + { + return Expansion(self(set), expansion, do_expansion); + } + +} // namespace samurai diff --git a/include/samurai/subset/lca_view.hpp b/include/samurai/subset/lca_view.hpp new file mode 100644 index 000000000..cf82b3de5 --- /dev/null +++ b/include/samurai/subset/lca_view.hpp @@ -0,0 +1,146 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_base.hpp" +#include "traversers/lca_traverser.hpp" + +namespace samurai +{ + + template + class LCAView; + + template + struct SetTraits> + { + static_assert(std::same_as, LCA>); + + template + using traverser_t = LCATraverser; + + static constexpr std::size_t dim() + { + return LCA::dim; + } + }; + + template + class LCAView : public SetBase> + { + using Self = LCAView; + + public: + + SAMURAI_SET_TYPEDEFS + + using const_interval_iterator = typename std::vector::const_iterator; + + explicit LCAView(const LCA& lca) + : m_lca(lca) + { + } + + inline std::size_t level_impl() const + { + return m_lca.level(); + } + + inline bool exist_impl() const + { + return !empty_impl(); + } + + inline bool empty_impl() const + { + return m_lca.empty(); + } + + template + inline void init_get_traverser_work_impl(const std::size_t, std::integral_constant) const + { + } + + template + inline void clear_get_traverser_work_impl(std::integral_constant) const + { + } + + template + inline traverser_t get_traverser_impl(const index_t& index, std::integral_constant d_ic) const + { + return get_traverser_impl_detail(index, + m_lca[Base::dim - 1].cbegin(), + m_lca[Base::dim - 1].cend(), + d_ic, + std::integral_constant{}); + } + + private: + + template + inline traverser_t get_traverser_impl_detail(const index_t& index, + const_interval_iterator begin_y_interval, + const_interval_iterator end_y_interval, + std::integral_constant d_ic, + std::integral_constant dCur_ic) const + { + if constexpr (dCur != Base::dim - 1) + { + const auto& y = index[dCur]; + const auto& y_offsets = m_lca.offsets(dCur + 1); + // we need to find an interval that contains y. + const auto y_interval_it = std::find_if(begin_y_interval, + end_y_interval, + [y](const auto& y_interval) + { + return y_interval.contains(y); + }); + if (y_interval_it != end_y_interval) + { + const std::size_t y_offset_idx = std::size_t(y + y_interval_it->index); + + const_interval_iterator begin_x_interval = m_lca[dCur].cbegin() + ptrdiff_t(y_offsets[y_offset_idx]); + const_interval_iterator end_x_interval = m_lca[dCur].cbegin() + ptrdiff_t(y_offsets[y_offset_idx + 1]); + + return get_traverser_impl_detail_wrapper(index, begin_x_interval, end_x_interval, d_ic, dCur_ic); + } + else + { + return get_traverser_impl_detail_wrapper(index, m_lca[dCur].cend(), m_lca[dCur].cend(), d_ic, dCur_ic); + } + } + else + { + return get_traverser_impl_detail_wrapper(index, m_lca[dCur].cbegin(), m_lca[dCur].cend(), d_ic, dCur_ic); + } + } + + template + inline traverser_t get_traverser_impl_detail_wrapper(const index_t& index, + const_interval_iterator begin_x_interval, + const_interval_iterator end_x_interval, + std::integral_constant d_ic, + std::integral_constant) const + { + if constexpr (d == dCur) + { + return traverser_t(begin_x_interval, end_x_interval); + } + else + { + return get_traverser_impl_detail(index, begin_x_interval, end_x_interval, d_ic, std::integral_constant{}); + } + } + + const LCA& m_lca; + }; + + template + LCAView> self(const LevelCellArray& lca) + { + return LCAView>(lca); + } + +} // namespace samurai diff --git a/include/samurai/subset/nary_set_operator.hpp b/include/samurai/subset/nary_set_operator.hpp new file mode 100644 index 000000000..2df60d937 --- /dev/null +++ b/include/samurai/subset/nary_set_operator.hpp @@ -0,0 +1,202 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "nary_set_operator_types.hpp" +#include "set_base.hpp" +#include "traversers/difference_id_traverser.hpp" +#include "traversers/difference_traverser.hpp" +#include "traversers/intersection_traverser.hpp" +#include "traversers/union_traverser.hpp" +#include "utils.hpp" + +namespace samurai +{ + enum class SetOperator + { + UNION, + INTERSECTION, + DIFFERENCE + }; + + template + class NArySetOperator; + + template + struct SetTraits> + { + static_assert((IsSet::value and ...)); + + template + using traverser_t = UnionTraverser...>; + + static constexpr std::size_t dim() + { + return std::tuple_element_t<0, std::tuple>::dim; + } + }; + + template + struct SetTraits> + { + static_assert((IsSet::value and ...)); + + template + using traverser_t = IntersectionTraverser...>; + + static constexpr std::size_t dim() + { + return std::tuple_element_t<0, std::tuple>::dim; + } + }; + + template + struct SetTraits> + { + static_assert((IsSet::value and ...)); + + template + using traverser_t = std::conditional_t...>, + DifferenceIdTraverser...>>; + + static constexpr std::size_t dim() + { + return std::tuple_element_t<0, std::tuple>::dim; + } + }; + + template + class NArySetOperator : public SetBase> + { + using Self = NArySetOperator; + + public: + + SAMURAI_SET_TYPEDEFS + + using Childrens = std::tuple; + + static constexpr std::size_t nIntervals = std::tuple_size_v; + + explicit NArySetOperator(const Sets&... sets) + : m_sets(sets...) + { + m_level = std::apply( + [](const auto&... set) -> std::size_t + { + return vmax(set.level()...); + }, + m_sets); + + enumerate_const_items(m_sets, + [this](const auto i, const auto& set) + { + m_shifts[i] = std::size_t(m_level - set.level()); + }); + } + + inline std::size_t level_impl() const + { + return m_level; + } + + inline bool exist_impl() const + { + return std::apply( + [](const auto first_set, const auto&... other_sets) -> std::size_t + { + if constexpr (op == SetOperator::UNION) + { + return first_set.exist() || (other_sets.exist() || ...); + } + else if constexpr (op == SetOperator::INTERSECTION) + { + return first_set.exist() && (other_sets.exist() && ...); + } + else + { + return first_set.exist(); + } + }, + m_sets); + } + + inline bool empty_impl() const + { + return Base::empty_default_impl(); + } + + template + inline void init_get_traverser_work_impl(const std::size_t n_traversers, std::integral_constant d_ic) const + { + static_for<0, nIntervals>::apply( + [this, n_traversers, d_ic](const auto i) -> void + { + std::get(m_sets).init_get_traverser_work(n_traversers, d_ic); + }); + } + + template + inline void clear_get_traverser_work_impl(std::integral_constant d_ic) const + { + static_for<0, nIntervals>::apply( + [this, d_ic](const auto i) -> void + { + std::get(m_sets).clear_get_traverser_work(d_ic); + }); + } + + template + inline traverser_t get_traverser_impl(const index_t& index, std::integral_constant d_ic) const + { + return get_traverser_impl_detail(index, d_ic, std::make_index_sequence{}); + } + + private: + + template + traverser_t + get_traverser_impl_detail(const index_t& index, std::integral_constant d_ic, std::index_sequence) const + { + return traverser_t(m_shifts, std::get(m_sets).get_traverser(index >> m_shifts[Is], d_ic)...); + } + + Childrens m_sets; + std::size_t m_level; + std::array m_shifts; + }; + + //////////////////////////////////////////////////////////////////////// + //// functions + //////////////////////////////////////////////////////////////////////// + + template + auto union_(const Sets&... sets) + requires(sizeof...(Sets) >= 1) + { + using Union = NArySetOperator...>; + + return Union(self(sets)...); + } + + template + auto intersection(const Sets&... sets) + requires(sizeof...(Sets) >= 2) + { + using Intersection = NArySetOperator...>; + + return Intersection(self(sets)...); + } + + template + auto difference(const Sets&... sets) + requires(sizeof...(Sets) >= 2) + { + using Difference = NArySetOperator...>; + + return Difference(self(sets)...); + } + +} // namespace samurai diff --git a/include/samurai/subset/node.hpp b/include/samurai/subset/node.hpp index 6114bf0f4..d36aa8009 100644 --- a/include/samurai/subset/node.hpp +++ b/include/samurai/subset/node.hpp @@ -3,687 +3,13 @@ #pragma once -#include -#include -#include -#include -#include -#include - -#include - -#include "../algorithm.hpp" #include "apply.hpp" -#include "concepts.hpp" -#include "samurai/list_of_intervals.hpp" -#include "start_end_fct.hpp" -#include "visitor.hpp" - -namespace samurai -{ - - template - void apply(Set&& global_set, Func&& func); - - template - class Subset - { - public: - - static constexpr std::size_t dim = get_set_dim_v; - using set_type = std::tuple; - using interval_t = get_interval_t; - - Subset(Op&& op, StartEndOp&& start_end_op, S&&... s) - : m_operator(std::forward(op)) - , m_start_end_op(std::forward(start_end_op)) - , m_s(std::forward(s)...) - , m_ref_level(compute_max(s.ref_level()...)) - , m_level(compute_max(s.level()...)) - , m_min_level(m_level) - { - std::apply( - [this](auto&&... args) - { - (args.ref_level(m_ref_level), ...); - }, - m_s); - } - - auto& on(auto level) - { - auto ilevel = static_cast(level); - if (ilevel > m_ref_level) - { - ref_level(ilevel); - } - m_min_level = std::min(m_min_level, ilevel); - m_level = ilevel; - return *this; - } - - template - void operator()(Func&& func) - { - apply(*this, std::forward(func)); - } - - template - void apply_op(ApplyOp&&... op) - { - auto func = [&](auto& interval, auto& index) - { - (op(m_level, interval, index), ...); - }; - apply(*this, func); - } - - inline void to_stream(std::ostream& os) - { - apply_op( - [&](const auto level, const auto& i, const auto& index) - { - os << "level: " << level << ", i: " << i << ", index: " << index << std::endl; - }); - } - - template - auto get_local_set(auto level, auto& index, Func_goback_beg&& goback_fct_beg, Func_goback_end&& goback_fct_end) - { - int shift = static_cast(this->ref_level()) - static_cast(this->level()); - m_start_end_op(m_level, m_min_level, m_ref_level); - - return std::apply( - [this, &index, shift, level, &goback_fct_beg, &goback_fct_end](auto&&... args) - { - return SetTraverser(shift, - get_operator(m_operator), - args.template get_local_set( - level, - index, - m_start_end_op.template goback(std::forward(goback_fct_beg)), - m_start_end_op.template goback(std::forward(goback_fct_end)))...); - }, - m_s); - } - - template - auto get_local_set(auto level, auto& index) - { - return get_local_set(level, index, default_function_(), default_function_()); - } - - template - auto get_start_and_stop_function(Func_start&& start_fct, Func_end&& end_fct) - { - m_start_end_op(m_level, m_min_level, m_ref_level); - - return std::apply( - [this, &start_fct, &end_fct](auto&& arg, auto&&... args) - { - if constexpr (std::is_same_v) - { - return std::make_tuple(std::move(arg.template get_start_and_stop_function( - m_start_end_op.template start(std::forward(start_fct)), - m_start_end_op.template end(std::forward(end_fct)))), - std::move(args.template get_start_and_stop_function( - m_start_end_op.template start(std::forward(start_fct)), - m_start_end_op.template end(std::forward(end_fct))))...); - } - else - { - return std::make_tuple(std::move(arg.template get_start_and_stop_function( - m_start_end_op.template start(std::forward(start_fct)), - m_start_end_op.template end(std::forward(end_fct)))), - std::move(args.template get_start_and_stop_function( - m_start_end_op.template start(std::forward(start_fct)), - m_start_end_op.template end(std::forward(end_fct))))...); - } - }, - m_s); - } - - template - auto get_start_and_stop_function() - { - return get_start_and_stop_function(default_function(), default_function()); - } - - auto level() const - { - return m_level; - } - - auto ref_level() const - { - return m_ref_level; - } - - void ref_level(auto level) - { - m_ref_level = level; - std::apply( - [this](auto&&... args) - { - (args.ref_level(m_ref_level), ...); - }, - m_s); - } - - bool empty() - { - return empty_check(*this); - } - - bool exist() const - { - return std::apply( - [this](auto&&... args) - { - return m_operator.exist(args...); - }, - m_s); - } - - protected: - - Op m_operator; - StartEndOp m_start_end_op; - set_type m_s; - std::size_t m_ref_level; - std::size_t m_level; - std::size_t m_min_level; - }; - - template - requires IsLCA - struct Self - { - static constexpr std::size_t dim = lca_t::dim; - using interval_t = typename lca_t::interval_t; - using value_t = typename interval_t::value_t; - - explicit Self(const lca_t& lca) - : m_lca(lca) - , m_level(lca.level()) - , m_ref_level(m_level) - , m_min_level(m_level) - { - } - - template - void operator()(Func&& func) - { - apply(*this, std::forward(func)); - } - - template - void apply_op(ApplyOp&&... op) - { - auto func = [&](auto& interval, auto& index) - { - (op(m_level, interval, index), ...); - }; - apply(*this, func); - } - - template - auto get_local_set(auto level, auto& index, Func_goback_beg&& goback_fct_beg, Func_goback_end&& goback_fct_end) - { - if (m_lca[d - 1].empty()) - { - return IntervalListVisitor(IntervalListRange(m_lca[d - 1], 0, 0)); - } - - if constexpr (dim == d) - { - m_offsets[d - 1].clear(); - m_offsets[d - 1].push_back({0, m_lca[d - 1].size()}); - - return IntervalListVisitor(m_lca.level(), - m_level, - m_ref_level, - IntervalListRange(m_lca[d - 1], 0, static_cast(m_lca[d - 1].size()))); - } - else - { - if (m_offsets[d].empty() || m_lca[d].empty()) - { - return IntervalListVisitor(IntervalListRange(m_lca[d - 1], 0, 0)); - } - - auto new_goback_fct_beg = m_func.template goback(std::forward(goback_fct_beg)); - - if (level <= m_level && level >= m_lca.level()) - { - m_offsets[d - 1].clear(); - - auto current_index = start_shift(new_goback_fct_beg(level, index[d - 1]).second, - static_cast(m_lca.level()) - static_cast(m_level)); - auto j = find_on_dim(m_lca, d, m_offsets[d][0][0], m_offsets[d][0][1], current_index); - - if (j == std::numeric_limits::max()) - { - return IntervalListVisitor(IntervalListRange(m_lca[d - 1], 0, 0)); - } - - auto io = static_cast(m_lca[d][j].index + current_index); - auto& offsets = m_lca.offsets(d); - m_offsets[d - 1].push_back({offsets[io], offsets[io + 1]}); - - return IntervalListVisitor(m_lca.level(), - m_level, - m_ref_level, - IntervalListRange(m_lca[d - 1], - static_cast(offsets[io]), - static_cast(offsets[io + 1]))); - } - else - { - auto new_goback_fct_end = m_func.template goback(std::forward(goback_fct_end)); - - auto min_index = start_shift(new_goback_fct_beg(level, index[d - 1]).second, - static_cast(m_lca.level()) - static_cast(m_level)); - - auto max_index = end_shift(new_goback_fct_end(level, index[d - 1] + 1).second, - static_cast(m_lca.level()) - static_cast(m_level)); - - m_work[d - 1].clear(); - m_offsets[d - 1].clear(); - - if constexpr (d == dim - 1) - { - auto j_min = lower_bound_interval(m_lca[d].begin() + static_cast(m_offsets[d][0][0]), - m_lca[d].begin() + static_cast(m_offsets[d][0][1]), - min_index); - auto j_max = upper_bound_interval(j_min, m_lca[d].begin() + static_cast(m_offsets[d][0][1]), max_index) - - 1; - - if (j_min != m_lca[d].end() && j_min <= j_max) - { - auto start_offset = static_cast(j_min->index + j_min->start); - if (j_min->contains(min_index)) - { - start_offset = static_cast(j_min->index + min_index); - } - - auto end_offset = static_cast(j_max->index + j_max->end); - if (j_max->contains(max_index)) - { - end_offset = static_cast(j_max->index + max_index); - } - - if (start_offset == end_offset) - { - return IntervalListVisitor(IntervalListRange(m_lca[d - 1], 0, 0)); - } - - m_offsets[d - 1].push_back({start_offset, end_offset}); - - ListOfIntervals list_of_intervals; - for (std::size_t o = m_lca.offsets(d)[start_offset]; o < m_lca.offsets(d)[end_offset]; ++o) - { - auto start = m_lca[d - 1][o].start; - auto end = m_lca[d - 1][o].end; - list_of_intervals.add_interval({start, end}); - } - - for (auto& i : list_of_intervals) - { - m_work[d - 1].push_back(i); - } - } - } - else - { - ListOfIntervals list_of_intervals; - - for (auto& offset : m_offsets[d]) - { - for (std::size_t io = offset[0]; io < offset[1]; ++io) - { - auto j_min = lower_bound_interval( - m_lca[d].begin() + static_cast(m_lca.offsets(d + 1)[io]), - m_lca[d].begin() + static_cast(m_lca.offsets(d + 1)[io + 1]), - min_index); - auto j_max = upper_bound_interval( - j_min, - m_lca[d].begin() + static_cast(m_lca.offsets(d + 1)[io + 1]), - max_index) - - 1; - - if (j_min != m_lca[d].begin() + static_cast(m_lca.offsets(d + 1)[io + 1]) && j_min <= j_max) - { - auto start_offset = static_cast(j_min->index + j_min->start); - if (j_min->contains(min_index)) - { - start_offset = static_cast(j_min->index + min_index); - } - - auto end_offset = static_cast(j_max->index + j_max->end); - if (j_max->contains(max_index)) - { - end_offset = static_cast(j_max->index + max_index); - } - - if (start_offset == end_offset) - { - continue; - } - - m_offsets[d - 1].push_back({start_offset, end_offset}); - - for (std::size_t o = m_lca.offsets(d)[start_offset]; o < m_lca.offsets(d)[end_offset]; ++o) - { - auto start = m_lca[d - 1][o].start; - auto end = m_lca[d - 1][o].end; - list_of_intervals.add_interval({start, end}); - } - } - } - - for (auto& i : list_of_intervals) - { - m_work[d - 1].push_back(i); - } - } - } - if (m_work[d - 1].empty()) - { - return IntervalListVisitor(IntervalListRange(m_lca[d - 1], 0, 0)); - } - return IntervalListVisitor(m_lca.level(), m_level, m_ref_level, IntervalListRange(m_lca[d - 1], m_work[d - 1])); - } - } - } - - template - auto get_local_set(auto level, auto& index) - - { - return get_local_set(level, index, default_function_(), default_function_()); - } - - template - auto get_start_and_stop_function(Func_start&& start_fct, Func_end&& end_fct) - { - m_func(m_level, m_min_level, m_ref_level); - auto new_start_fct = m_func.template start(std::forward(start_fct)); - auto new_end_fct = m_func.template end(std::forward(end_fct)); - return std::make_tuple(std::move(new_start_fct), std::move(new_end_fct)); - } - - template - auto get_start_and_stop_function() - - { - return get_start_and_stop_function(default_function(), default_function()); - } - - auto ref_level() const - { - return m_ref_level; - } - - void ref_level(auto level) - { - m_ref_level = level; - } - - auto level() const - { - return m_level; - } - - auto& on(auto level) - { - m_ref_level = std::max(m_ref_level, static_cast(level)); - m_min_level = std::min(m_min_level, static_cast(level)); - m_level = static_cast(level); - return *this; - } - - bool exist() const - { - return !m_lca.empty(); - } - - bool empty() const - { - return !m_lca.empty(); - } - - const lca_t& m_lca; - std::size_t m_level; - std::size_t m_ref_level; - std::size_t m_min_level; - start_end_function m_func; - std::array, dim - 1> m_work; - std::array>, dim> m_offsets; - }; - - namespace detail - { - template - auto transform(const LevelCellArray& lca) - { - return Self>(lca); - } - - template - auto transform(LevelCellArray& lca) - { - return Self>(lca); - } - - template - auto transform(E&& e) - { - return std::forward(e); - } - } - - template - inline std::ostream& operator<<(std::ostream& out, Subset& subset) - { - subset.to_stream(out); - return out; - } - - template - auto intersection(sets_t&&... sets) - { - static constexpr std::size_t dim = get_set_dim_v; - return std::apply( - [](auto&&... args) - { - return Subset(IntersectionOp(), start_end_function(), std::forward(args)...); - }, - std::make_tuple(detail::transform(std::forward(sets))...)); - } - - template - auto union_(sets_t&&... sets) - { - static constexpr std::size_t dim = get_set_dim_v; - return std::apply( - [](auto&&... args) - { - return Subset(UnionOp(), start_end_function(), std::forward(args)...); - }, - std::make_tuple(detail::transform(std::forward(sets))...)); - } - - template - auto difference(sets_t&&... sets) - { - static constexpr std::size_t dim = get_set_dim_v; - return std::apply( - [](auto&&... args) - { - return Subset(DifferenceOp(), start_end_function(), std::forward(args)...); - }, - std::make_tuple(detail::transform(std::forward(sets))...)); - } - - template - auto translate(set_t&& set, const stencil_t& stencil) - { - constexpr std::size_t dim = std::decay_t::dim; - return Subset(SelfOp(), - start_end_translate_function(xt::xtensor_fixed>(stencil)), - detail::transform(std::forward(set))); - } - - template - auto contraction(set_t&& set, int c) - { - constexpr std::size_t dim = std::decay_t::dim; - return Subset(SelfOp(), start_end_contraction_function(c), detail::transform(std::forward(set))); - } - - template - auto self(lca_t&& lca) - { - return Self>(std::forward(lca)); - } - - //----------------------------------------------------------------// - // Contract // - //----------------------------------------------------------------// - - template - auto contract_rec(const SubsetOrLCA& set, int width, const std::array& contract_directions) - { - static constexpr std::size_t dim = SubsetOrLCA::dim; - if constexpr (direction_index < dim) - { - using direction_t = xt::xtensor_fixed>; - - auto contracted_in_other_dirs = contract_rec(set, width, contract_directions); - direction_t dir; - dir.fill(0); - dir[direction_index] = contract_directions[direction_index] ? width : 0; - - return intersection(contracted_in_other_dirs, translate(set, dir), translate(set, -dir)); - } - else - { - if constexpr (IsLCA) - { - return self(set); - } - else - { - return set; - } - } - } - - /** - * @brief Contract a set in the specified directions. - * - * @tparam SubsetOrLCA The type of the set to contract. - * @param set The set or LevelCellArray to contract. - * @param width The contraction width. - * @param contract_directions An array indicating which directions to contract (true for contraction, false for no contraction). - * @return A new set that is contracted in the specified directions. - */ - template - auto contract(const SubsetOrLCA& set, std::size_t width, const std::array& contract_directions) - { - return contract_rec<0>(set, static_cast(width), contract_directions); - } - - /** - * @brief Contract a set in all directions. - * - * This function is a convenience wrapper that contracts the set in all dimensions. - * - * @tparam SubsetOrLCA The type of the set to contract. - * @param set The set or LevelCellArray to contract. - * @param width The contraction width. - * @return A new set that is contracted in all directions. - */ - template - auto contract(const SubsetOrLCA& set, std::size_t width) - { - std::array contract_directions; - std::fill(contract_directions.begin(), contract_directions.end(), true); - return contract(set, width, contract_directions); - } - - //--------------------------------------------------------------// - // Expand // - //--------------------------------------------------------------// - - template - auto expand_rec(const SubsetOrLCA& set, const std::array& expand_directions) - { - static constexpr std::size_t dim = SubsetOrLCA::dim; - using direction_t = xt::xtensor_fixed>; - - if constexpr (width == 0 && direction_index == dim) - { - if constexpr (IsLCA) - { - return self(set); - } - else - { - return set; - } - } - else if constexpr (width > 0) - { - auto expanded_layer = expand_rec(set, expand_directions); - direction_t dir; - dir.fill(0); - dir[direction_index] = expand_directions[direction_index] ? width : 0; - - return union_(expanded_layer, translate(set, dir), translate(set, -dir)); - } - else if constexpr (direction_index < dim) - { - auto expanded_in_other_dirs = expand_rec(set, expand_directions); - direction_t dir; - dir.fill(0); - dir[direction_index] = expand_directions[direction_index] ? width : 0; - - return union_(expanded_in_other_dirs, translate(set, dir), translate(set, -dir)); - } - } - - /** - * @brief Expand a set in the specified directions. - * - * @tparam SubsetOrLCA The type of the set to expand. - * @param set The set or LevelCellArray to expand. - * @param width The expansion width. - * @param expand_directions An array indicating which directions to expand (true for expansion, false for no expansion). - * @return A new set that is expanded in the specified directions. - */ - template - auto expand(const SubsetOrLCA& set, const std::array& expand_directions) - { - return expand_rec<0, width>(set, expand_directions); - } - - /** - * @brief Expand a set in all Cartesian directions (no diagonals!). - * - * This function is a convenience wrapper that expands the set in all dimensions. - * - * @tparam SubsetOrLCA The type of the set to expand. - * @param set The set or LevelCellArray to expand. - * @param width The expansion width. - * @return A new set that is expanded in all directions. - */ - template - auto expand(const SubsetOrLCA& set) - { - std::array expand_directions; - std::fill(expand_directions.begin(), expand_directions.end(), true); - return expand(set, expand_directions); - } -} +#include "box_view.hpp" +#include "contraction.hpp" +#include "expansion.hpp" +#include "lca_view.hpp" +#include "nary_set_operator.hpp" +#include "projection.hpp" +#include "set_base.hpp" +#include "translation.hpp" +#include "utils.hpp" diff --git a/include/samurai/subset/projection.hpp b/include/samurai/subset/projection.hpp new file mode 100644 index 000000000..d651e25e2 --- /dev/null +++ b/include/samurai/subset/projection.hpp @@ -0,0 +1,187 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../samurai_config.hpp" +#include "../static_algorithm.hpp" +#include "set_base.hpp" +#include "traversers/projection_traverser.hpp" + +namespace samurai +{ + + template + class Projection; + + template + struct SetTraits> + { + static_assert(IsSet::value); + + template + using traverser_t = ProjectionTraverser>; + + static constexpr std::size_t dim() + { + return Set::dim; + } + }; + + namespace detail + { + template + struct ProjectionWork; + + template + struct ProjectionWork> + { + template + using child_traverser_t = typename Set::template traverser_t; + + template + using array_of_child_traverser_t = std::vector>; + + using Type = std::tuple...>; + }; + } // namespace detail + + template + class Projection : public SetBase> + { + using Self = Projection; + using ChildTraverserArray = detail::ProjectionWork>::Type; + + public: + + SAMURAI_SET_TYPEDEFS + + Projection(const Set& set, const std::size_t level) + : m_set(set) + , m_level(level) + { + if (m_level < m_set.level()) + { + m_projectionType = ProjectionType::COARSEN; + m_shift = m_set.level() - m_level; + } + else + { + m_projectionType = ProjectionType::REFINE; + m_shift = m_level - m_set.level(); + } + } + + // we need to define a custom copy and move constructor because + // we do not want to copy m_work_offsetRanges + Projection(const Projection& other) + : m_set(other.m_set) + , m_level(other.m_level) + , m_projectionType(other.m_projectionType) + , m_shift(other.m_shift) + { + } + + Projection(Projection&& other) + : m_set(std::move(other.m_set)) + , m_level(std::move(other.m_level)) + , m_projectionType(std::move(other.m_projectionType)) + , m_shift(std::move(other.m_shift)) + { + } + + inline std::size_t level_impl() const + { + return m_level; + } + + inline bool exist_impl() const + { + return m_set.exist(); + } + + inline bool empty_impl() const + { + return m_set.empty(); + } + + template + inline void init_get_traverser_work_impl(const std::size_t n_traversers, std::integral_constant d_ic) const + { + const std::size_t my_work_size_per_traverser = (m_projectionType == ProjectionType::COARSEN and d != Base::dim - 1) + ? (1 << m_shift) + : 1; + const std::size_t my_work_size = n_traversers * my_work_size_per_traverser; + + auto& childTraversers = std::get(m_work_childTraversers); + childTraversers.clear(); + childTraversers.reserve(my_work_size); + + m_set.init_get_traverser_work(my_work_size, d_ic); + } + + template + inline void clear_get_traverser_work_impl(std::integral_constant d_ic) const + { + auto& childTraversers = std::get(m_work_childTraversers); + + childTraversers.clear(); + + m_set.clear_get_traverser_work(d_ic); + } + + template + inline traverser_t get_traverser_impl(const index_t& _index, std::integral_constant d_ic) const + { + auto& childTraversers = std::get(m_work_childTraversers); + + if (m_projectionType == ProjectionType::COARSEN) + { + if constexpr (d != Base::dim - 1) + { + const std::size_t old_capacity = childTraversers.capacity(); + + const value_t ymin = _index[d] << m_shift; + const value_t ybound = (_index[d] + 1) << m_shift; + + xt::xtensor_fixed> index(_index << m_shift); + + const auto childTraversers_begin = childTraversers.end(); + + for (index[d] = ymin; index[d] != ybound; ++index[d]) + { + childTraversers.push_back(m_set.get_traverser(index, d_ic)); + if (childTraversers.back().is_empty()) + { + childTraversers.pop_back(); + } + } + + assert(childTraversers.capacity() == old_capacity); + + return traverser_t(childTraversers_begin, childTraversers.end(), m_shift); + } + else + { + childTraversers.push_back(m_set.get_traverser(_index << m_shift, d_ic)); + return traverser_t(std::prev(childTraversers.end()), m_projectionType, m_shift); + } + } + else + { + childTraversers.push_back(m_set.get_traverser(_index >> m_shift, d_ic)); + return traverser_t(std::prev(childTraversers.end()), m_projectionType, m_shift); + } + } + + private: + + Set m_set; + std::size_t m_level; + ProjectionType m_projectionType; + std::size_t m_shift; + + mutable ChildTraverserArray m_work_childTraversers; + }; + +} // namespace samurai diff --git a/include/samurai/subset/projection_loi.hpp b/include/samurai/subset/projection_loi.hpp new file mode 100644 index 000000000..2e13c6676 --- /dev/null +++ b/include/samurai/subset/projection_loi.hpp @@ -0,0 +1,202 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../samurai_config.hpp" +#include "../static_algorithm.hpp" +#include "set_base.hpp" +#include "traversers/last_dim_projection_loi_traverser.hpp" +#include "traversers/projection_loi_traverser.hpp" + +#include + +namespace samurai +{ + + template + class ProjectionLOI; + + template + struct SetTraits> + { + static_assert(IsSet::value); + + template + using child_traverser_t = typename Set::template traverser_t; + + template + using traverser_t = std::conditional_t>, + ProjectionLOITraverser>>; + + static constexpr std::size_t dim() + { + return Set::dim; + } + }; + + namespace detail + { + template + struct ProjectionLOIWork; + + template + struct ProjectionLOIWork> + { + template + using child_traverser_t = typename Set::template traverser_t; + + template + using work_t = ListOfIntervals::value_t>; + + using Type = std::tuple...>; + }; + } // namespace detail + + template + class ProjectionLOI : public SetBase> + { + using Self = ProjectionLOI; + using ListOfIntervals = detail::ProjectionLOIWork>::Type; + + public: + + SAMURAI_SET_TYPEDEFS + + using Index = xt::xtensor_fixed>; + + ProjectionLOI(const Set& set, const std::size_t level) + : m_set(set) + , m_level(level) + { + if (m_level < m_set.level()) + { + m_projectionType = ProjectionType::COARSEN; + m_shift = m_set.level() - m_level; + } + else + { + m_projectionType = ProjectionType::REFINE; + m_shift = m_level - m_set.level(); + } + } + + // we need to define a custom copy and move constructor because + // we do not want to copy m_work_offsetRanges + ProjectionLOI(const ProjectionLOI& other) + : m_set(other.m_set) + , m_level(other.m_level) + , m_projectionType(other.m_projectionType) + , m_shift(other.m_shift) + { + } + + ProjectionLOI(ProjectionLOI&& other) + : m_set(std::move(other.m_set)) + , m_level(std::move(other.m_level)) + , m_projectionType(std::move(other.m_projectionType)) + , m_shift(std::move(other.m_shift)) + { + } + + inline std::size_t level_impl() const + { + return m_level; + } + + inline bool exist_impl() const + { + return m_set.exist(); + } + + inline bool empty_impl() const + { + return m_set.empty(); + } + + template + inline void init_get_traverser_work_impl(const std::size_t n_traversers, std::integral_constant d_ic) const + { + assert(n_traversers == 1); + + m_set.init_get_traverser_work(n_traversers, d_ic); + } + + template + inline void clear_get_traverser_work_impl(std::integral_constant d_ic) const + { + m_set.clear_get_traverser_work(d_ic); + } + + template + inline traverser_t get_traverser_impl(const index_t& index, std::integral_constant d_ic) const + { + auto& listOfIntervals = std::get(m_work_listOfIntervals); + listOfIntervals.clear(); + + if (m_projectionType == ProjectionType::COARSEN) + { + if constexpr (d != Base::dim - 1) + { + Index index_min(index << m_shift); + Index index_max((index + 1) << m_shift); + + Index index_rec; + fill_list_of_interval_rec(index_min, index_max, index_rec, d_ic, std::integral_constant{}); + + return traverser_t(m_set.get_traverser(index << m_shift, d_ic), listOfIntervals.cbegin(), listOfIntervals.cend()); + } + else + { + return traverser_t(m_set.get_traverser(index << m_shift, d_ic), m_projectionType, m_shift); + } + } + else + { + return traverser_t(m_set.get_traverser(index >> m_shift, d_ic), m_projectionType, m_shift); + } + } + + private: + + template + inline void fill_list_of_interval_rec(const Index& index_min, + const Index& index_max, + index_t& index, + std::integral_constant d_ic, + std::integral_constant dCur_ic) const + { + using child_traverser_t = typename Set::template traverser_t; + using child_current_interval_t = typename child_traverser_t::current_interval_t; + + for (child_traverser_t traverser = m_set.get_traverser(index, dCur_ic); !traverser.is_empty(); traverser.next_interval()) + { + child_current_interval_t interval = traverser.current_interval(); + + if constexpr (dCur == d) + { + std::get(m_work_listOfIntervals).add_interval(interval >> m_shift); + } + else + { + const auto index_start = std::max(interval.start, index_min[dCur - 1]); + const auto index_end = std::min(interval.end, index_max[dCur - 1]); + + for (index[dCur - 1] = index_start; index[dCur - 1] < index_end; ++index[dCur - 1]) + { + fill_list_of_interval_rec(index_min, index_max, index, d_ic, std::integral_constant{}); + } + } + } + } + + Set m_set; + std::size_t m_level; + ProjectionType m_projectionType; + std::size_t m_shift; + + mutable ListOfIntervals m_work_listOfIntervals; + }; + +} // namespace samurai diff --git a/include/samurai/subset/projection_type.hpp b/include/samurai/subset/projection_type.hpp new file mode 100644 index 000000000..5489a0c0b --- /dev/null +++ b/include/samurai/subset/projection_type.hpp @@ -0,0 +1,13 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +namespace samurai +{ + enum class ProjectionType + { + COARSEN, + REFINE + }; +} // namespace samurai diff --git a/include/samurai/subset/range_traversers/box_traverser_range.hpp b/include/samurai/subset/range_traversers/box_traverser_range.hpp new file mode 100644 index 000000000..a7b5817fe --- /dev/null +++ b/include/samurai/subset/range_traversers/box_traverser_range.hpp @@ -0,0 +1,110 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../../interval.hpp" +#include "../traversers/box_traverser.hpp" +#include "set_traverser_range_base.hpp" + +namespace samurai +{ + template + class Box; + + template + class BoxTraverserRange; + + template + struct SetTraverserRangeTraits> + { + static_assert(std::same_as, B>); + + class Iterator + { + public: + + using index_t = typename interval_t::value_t; + + using iterator_category = std::forward_iterator_tag; + using difference_type = index_t; + using value_type = BoxTraverser; + using reference = BoxTraverser; + + Iterator(const interval_t& interval, const index_t& index) + : m_interval(interval) + , m_index(first_index) + { + } + + reference operator*() const + { + return reference(interval); + } + + Iterator& operator++() + { + ++m_index; + return *this; + } + + Iterator operator++(int) + { + Iterator tmp = *this; + ++(*this); + return tmp; + } + + friend bool operator==(const Iterator& a, const Iterator& b) + { + return a.m_index == b.m_index; + }; + + friend bool operator!=(const Iterator& a, const Iterator& b) + { + return a.m_index != b.m_index; + }; + + private: + + interval_t m_interval; + index_t m_index; + }; + }; + + template + class BoxTraverserRange : public SetTraverserRangeBase> + { + using Self = BoxTraverserRange; + + public: + + SAMURAI_SET_TRAVERSER_RANGE_TYPEDEFS + + using interval_t = Interval; + using index_t = typename interval_t::value_t; + + BoxTraverserRange(const interval_t& interval, const index_t& first_index, const index_t& last_index) + : m_interval(interval) + , m_first_index(first_index) + , m_last_index(last_index) + { + } + + Iterator begin_impl() + { + return Iterator(m_interval, m_first_index); + } + + Iterator end_impl() + { + return Iterator(m_interval, m_last_index); + } + + private: + + interval_t m_interval; + index_t m_first_index; + index_t m_last_index; + }; +} // namespace samurai diff --git a/include/samurai/subset/range_traversers/lca_traverser_range.hpp b/include/samurai/subset/range_traversers/lca_traverser_range.hpp new file mode 100644 index 000000000..e1fae88a9 --- /dev/null +++ b/include/samurai/subset/range_traversers/lca_traverser_range.hpp @@ -0,0 +1,114 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include +#include + +#include "lca_traverser_range_item.hpp" +#include "set_traverser_range_base.hpp" + +namespace samurai +{ + + template + class LevelCellArray; + + template + class LCATraverserRange; + + template + struct SetTraverserRangeTraits> + { + static_assert(std::same_as, LCA>); + + using interval_iterator = typename std::vector::const_iterator; + using offset_iterator = typename std::vector::iterator; + + class Iterator + { + public: + + using iterator_category = std::forward_iterator_tag; + using difference_type = std::ptrdiff_t; + using value_type = LCATraverserRangeItem; + using reference = LCATraverserRangeItem; + + Iterator(const interval_iterator first_interval, const offset_iterator offset) + : m_first_interval(first_interval) + , m_offset(offset) + { + } + + reference operator*() const + { + return reference(m_first_interval, m_offset); + } + + Iterator& operator++() + { + ++m_offset; + return *this; + } + + Iterator operator++(int) + { + Iterator tmp = *this; + ++(*this); + return tmp; + } + + friend bool operator==(const Iterator& a, const Iterator& b) + { + return a.m_first_interval == b.m_first_interval and a.m_offset == b.m_offset; + }; + + friend bool operator!=(const Iterator& a, const Iterator& b) + { + return a.m_ptr != b.m_ptr or a.m_offset != b.m_offset; + }; + + private: + + interval_iterator m_first_interval; + offset_iterator m_offset; + }; + }; + + template + class LCATraverserRange : public SetTraverserRangeBase> + { + using Self = LCATraverserRange; + + public: + + SAMURAI_SET_TRAVERSER_RANGE_TYPEDEFS + using interval_iterator = typename std::vector::const_iterator; + using offset_iterator = typename std::vector::iterator; + + LCATraverserRange(const interval_iterator first_interval, const offset_iterator first_offsets, const offset_iterator last_offsets) + : m_first_interval(first_interval) + , m_first_offsets(first_offsets) + , m_last_offsets(last_offsets) + { + } + + Iterator begin_impl() + { + return Iterator(m_first_interval, m_first_offsets); + } + + Iterator end_impl() + { + return Iterator(m_first_interval, std::prev(m_last_offsets)); + } + + private: + + interval_iterator m_first_interval; + offset_iterator m_first_offsets; + offset_iterator m_last_offsets; + }; + +} // namespace samurai diff --git a/include/samurai/subset/range_traversers/lca_traverser_range_item.hpp b/include/samurai/subset/range_traversers/lca_traverser_range_item.hpp new file mode 100644 index 000000000..99d6d855f --- /dev/null +++ b/include/samurai/subset/range_traversers/lca_traverser_range_item.hpp @@ -0,0 +1,67 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include +#include + +#include "../traversers/set_traverser_base.hpp" + +namespace samurai +{ + + template + class LevelCellArray; + + template + struct SetTraverserTraits> + { + static_assert(std::same_as, LCA>); + + using interval_t = typename LCA::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class LCATraverserRangeItem : public SetTraverserBase> + { + using Self = LCATraverserRangeItem; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + using interval_iterator = typename std::vector::const_iterator; + using offset_iterator = typename std::vector::iterator; + + LCATraverserRangeItem(const interval_iterator first_interval, const offset_iterator offset) + : m_first_interval(first_interval) + , m_offset(*offset) + , m_offset_bound(*(offset + 1)) + { + } + + inline bool is_empty_impl() const + { + return m_offset == m_offset_bound; + } + + inline void next_interval_impl() + { + ++m_offset; + } + + inline current_interval_t current_interval_impl() const + { + return *(m_first_interval + m_offset); + } + + private: + + interval_iterator m_first_interval; + std::size_t& m_offset; + std::size_t m_offset_bound; + }; + +} // namespace samurai diff --git a/include/samurai/subset/range_traversers/nary_traverser_range.hpp b/include/samurai/subset/range_traversers/nary_traverser_range.hpp new file mode 100644 index 000000000..e5345f7e7 --- /dev/null +++ b/include/samurai/subset/range_traversers/nary_traverser_range.hpp @@ -0,0 +1,120 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "nary_traverser_type.hpp" + +namespace samurai +{ + template + class NaryOperatorTraverserRange; + + template + struct SetTraverserRangeTraits> + { + static_assert((IsSetTraverserRange::value and ...)); + + class Iterator + { + public: + + using ChildrenIterators = std::tuple; + + using iterator_category = std::forward_iterator_tag; + using difference_type = std::ptrdiff_t; + using value_type = typename NAryTraverserType::Type; + using reference = value_type; + + static constexpr std::size_t nIntervals = std::tuple_size::value; + + Iterator(const ChildrenIterators& innerIterators) + : m_childrenIterators(innerIterators) + { + } + + reference operator*() const + { + return std::apply( + [](const auto&... childrenIterators) -> void + { + return reference((*childrenIterators)...); + }, + m_childrenIterators); + } + + Iterator& operator++() + { + static_for<0, nIntervals>::apply( + [this](const auto i) -> void + { + ++std::get(m_childrenIterators); + }); + return *this; + } + + Iterator operator++(int) + { + Iterator tmp = *this; + ++(*this); + return tmp; + } + + friend bool operator==(const Iterator& a, const Iterator& b) + { + return a.m_childrenIterators == b.m_childrenIterators; + }; + + friend bool operator!=(const Iterator& a, const Iterator& b) + { + return a.m_childrenIterators != b.m_childrenIterators; + }; + + private: + + ChildrenIterators m_childrenIterators; + }; + }; + + template + class NaryOperatorTraverserRange : public SetTraverserRangeBase> + { + using Self = TranslationTraverserRange; + + public: + + SAMURAI_SET_TRAVERSER_RANGE_TYPEDEFS + + using Childrens = std::tuple; + + TranslationTraverserRange(const SetTraverserRanges&... set_traverser_ranges) + : m_set_traverser_ranges(set_traverser_ranges...) + { + } + + Iterator begin_impl() + { + return std::apply( + [](auto&... innerIterators) -> Iterator + { + return Iterator(std::make_tuple(innerIterators.begin()...)); + }, + m_set_traverser_ranges); + } + + Iterator end_impl() + { + return std::apply( + [](auto&... innerIterators) -> Iterator + { + return Iterator(std::make_tuple(innerIterators.end()...)); + }, + m_set_traverser_ranges); + } + + private: + + Childrens m_set_traverser_ranges; + }; + +} // namespace samurai diff --git a/include/samurai/subset/range_traversers/nary_traverser_type.hpp b/include/samurai/subset/range_traversers/nary_traverser_type.hpp new file mode 100644 index 000000000..33ba899c0 --- /dev/null +++ b/include/samurai/subset/range_traversers/nary_traverser_type.hpp @@ -0,0 +1,48 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../difference_id_traverser.hpp" +#include "../difference_traverser.hpp" +#include "../intersection_traverser.hpp" +#include "../union_traverser.hpp" + +namespace samurai +{ + enum class NAryTraverserType + { + UNION, + INTERSECTION, + DIFFERENCE, + DIFFERENCE_ID + }; + + template + struct NAryTraverserTypeTraits; + + template + struct NAryTraverserTypeTraits + { + using Type = UnionTraverser; + }; + + template + struct NAryTraverserTypeTraits + { + using Type = IntersectionTraverser; + }; + + template + struct NAryTraverserTypeTraits + { + using Type = DifferenceTraverser; + }; + + template + struct NAryTraverserTypeTraits + { + using Type = DifferenceIdTraverser; + }; + +} // namespace samurai diff --git a/include/samurai/subset/range_traversers/set_traverser_range_base.hpp b/include/samurai/subset/range_traversers/set_traverser_range_base.hpp new file mode 100644 index 000000000..c741e0f3a --- /dev/null +++ b/include/samurai/subset/range_traversers/set_traverser_range_base.hpp @@ -0,0 +1,57 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include +#include + +namespace samurai +{ + template + struct SetTraverserRangeTraits; + + /* + * For the sake of conscision, the ranges are only forward range + * but they could easily be converted to a random_access range. + */ + template + class SetTraverserRangeBase + { + public: + + using DerivedTraits = SetTraverserRangeTraits; + + using Iterator = typename DerivedTraits::Iterator; + + const Derived& derived_cast() const + { + return static_cast(*this); + } + + Derived& derived_cast() + { + return static_cast(*this); + } + + Iterator begin() + { + derived_cast().begin_impl(); + } + + Iterator end() + { + derived_cast().end_impl(); + } + }; + + template + struct IsSetTraverserRange : std::bool_constant, T>::value> + { + }; + +#define SAMURAI_SET_TRAVERSER_RANGE_TYPEDEFS \ + using Base = SetTraverserBase; \ + using Iterator = typename Base::Iterator; + +} // namespace samurai diff --git a/include/samurai/subset/range_traversers/translation_traverser_range.hpp b/include/samurai/subset/range_traversers/translation_traverser_range.hpp new file mode 100644 index 000000000..b41f0a5b4 --- /dev/null +++ b/include/samurai/subset/range_traversers/translation_traverser_range.hpp @@ -0,0 +1,107 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_range_base.hpp" + +namespace samurai +{ + + template + class TranslationTraverserRange; + + template + struct SetTraverserRangeTraits> + { + static_assert(IsSetTraverserRange::value); + + using ChildIterator = typename SetTraverserRange::Iterator; + using const_ChildIterator = typename SetTraverserRange::constIterator; + + using Translation = typename ChildIterator::value_type::value_type; + + class Iterator + { + public: + + using iterator_category = std::forward_iterator_tag; + using difference_type = std::ptrdiff_t; + using value_type = TranslationTraverser; + using reference = value_type; + + Iterator(const InnerIterator innerIterator, const Translation& translation) + : m_innerIterator(innerIterator) + , m_translation(translation) + { + } + + reference operator*() const + { + return reference(*m_innerIterator, m_translation); + } + + Iterator& operator++() + { + ++m_innerIterator; + return *this; + } + + Iterator operator++(int) + { + Iterator tmp = *this; + ++(*this); + return tmp; + } + + friend bool operator==(const Iterator& a, const Iterator& b) + { + return a.m_innerIterator == b.m_innerIterator; + }; + + friend bool operator!=(const Iterator& a, const Iterator& b) + { + return a.m_innerIterator != b.m_innerIterator; + }; + + private: + + ChildIterator m_innerIterator; + Translation m_translation; + }; + }; + + template + class TranslationTraverserRange : public SetTraverserRangeBase> + { + using Self = TranslationTraverserRange; + + public: + + SAMURAI_SET_TRAVERSER_RANGE_TYPEDEFS + + using Translation = typename SetTraverserRangeTraits::Translation; + + TranslationTraverserRange(const SetTraverserRange& set_traverser_range, const Translation& translation) + : m_set_traverser_range(set_traverser_range) + , m_translation(translation) + { + } + + Iterator begin_impl() + { + return Iterator(m_set_traverser_range.begin(), m_translation); + } + + Iterator end_impl() + { + return Iterator(m_set_traverser_range.end(), m_translation); + } + + private: + + SetTraverserRange m_set_traverser_range; + Translation m_translation; + }; + +} // namespace samurai diff --git a/include/samurai/subset/set_base.hpp b/include/samurai/subset/set_base.hpp new file mode 100644 index 000000000..1f0aff9ca --- /dev/null +++ b/include/samurai/subset/set_base.hpp @@ -0,0 +1,205 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include + +#include + +#include "../samurai_config.hpp" +#include "traversers/set_traverser_base.hpp" + +namespace samurai +{ + //////////////////////////////////////////////////////////////////////// + //// Forward Declarations + //////////////////////////////////////////////////////////////////////// + + template + struct SetTraits; + + template + class SetBase; + + template + class Projection; + + template + class ProjectionLOI; + + template + void apply(const SetBase& set, Func&& func); + + //////////////////////////////////////////////////////////////////////// + //// Class definition + //////////////////////////////////////////////////////////////////////// + + template + class SetBase + { + using DerivedTraits = SetTraits; + + public: + + template + using traverser_t = typename DerivedTraits::template traverser_t; + using interval_t = typename traverser_t<0>::interval_t; + using value_t = typename interval_t::value_t; + + using to_lca_t = LevelCellArray; + using to_lca_coord_t = typename to_lca_t::coords_t; + + using ProjectionMethod = std::conditional_t, Projection>; + + static constexpr std::size_t dim = DerivedTraits::dim(); + + const Derived& derived_cast() const + { + return static_cast(*this); + } + + Derived& derived_cast() + { + return static_cast(*this); + } + + inline std::size_t level() const + { + return derived_cast().level_impl(); + } + + inline bool exist() const + { + return derived_cast().exist_impl(); + } + + inline bool empty() const + { + return derived_cast().empty_impl(); + } + + template + inline void init_get_traverser_work(const std::size_t n_traversers, std::integral_constant d_ic) const + { + derived_cast().init_get_traverser_work_impl(n_traversers, d_ic); + } + + template + inline void clear_get_traverser_work(std::integral_constant d_ic) const + { + derived_cast().clear_get_traverser_work_impl(d_ic); + } + + template + inline traverser_t get_traverser(const index_t& index, std::integral_constant d_ic) const + { + return derived_cast().get_traverser_impl(index, d_ic); + } + + inline ProjectionMethod on(const std::size_t level); + + template + void operator()(Func&& func) const + { + apply(derived_cast(), std::forward(func)); + } + + template + void apply_op(ApplyOp&&... op) const + { + const std::size_t l = level(); + + auto func = [l, &op...](auto& interval, auto& index) + { + (op(l, interval, index), ...); + }; + apply(derived_cast(), func); + } + + to_lca_t to_lca() const + { + return to_lca_t(*this); + } + + to_lca_t to_lca(const to_lca_coord_t& origin_point, const double scaling_factor) const + { + return to_lca_t(*this, origin_point, scaling_factor); + } + + protected: + + inline bool empty_default_impl() const + { + xt::xtensor_fixed> index; + return empty_default_impl_rec(index, std::integral_constant{}); + } + + template + bool empty_default_impl_rec(index_t& index, std::integral_constant d_ic) const + { + using current_interval_t = typename traverser_t::current_interval_t; + + init_get_traverser_work(1, d_ic); + + for (traverser_t traverser = get_traverser(index, d_ic); !traverser.is_empty(); traverser.next_interval()) + { + current_interval_t interval = traverser.current_interval(); + + if constexpr (d == 0) + { + return false; + } + else + { + for (index[d - 1] = interval.start; index[d - 1] != interval.end; ++index[d - 1]) + { + if (not empty_default_impl_rec(index, std::integral_constant{})) + { + return false; + } + } + } + } + + clear_get_traverser_work(d_ic); + + return true; + } + }; + +#define SAMURAI_SET_TYPEDEFS \ + using Base = SetBase; \ + \ + template \ + using traverser_t = typename Base::template traverser_t; \ + \ + using interval_t = typename Base::interval_t; \ + using value_t = typename Base::value_t; + + template + struct IsSet : std::bool_constant, T>::value> + { + }; + + template + const Set& self(const SetBase& set) + { + return set.derived_cast(); + } + +} // namespace samurai + +#include "projection.hpp" +#include "projection_loi.hpp" + +namespace samurai +{ + + template + auto SetBase::on(const std::size_t level) -> ProjectionMethod + { + return ProjectionMethod(derived_cast(), level); + } + +} diff --git a/include/samurai/subset/start_end_fct.hpp b/include/samurai/subset/start_end_fct.hpp deleted file mode 100644 index 6c304ee28..000000000 --- a/include/samurai/subset/start_end_fct.hpp +++ /dev/null @@ -1,275 +0,0 @@ -// Copyright 2018-2025 the samurai's authors -// SPDX-License-Identifier: BSD-3-Clause - -#pragma once - -#include "utils.hpp" - -namespace samurai -{ - - inline auto default_function() - { - return [](auto, auto i, auto) - { - return i; - }; - } - - inline auto default_function_() - { - return [](auto level, auto i) - { - return std::make_pair(level, i); - }; - } - - template - struct start_end_function - { - auto& operator()(std::size_t level, std::size_t min_level, std::size_t max_level) - { - m_level = level; - m_min_level = min_level; - m_shift = static_cast(max_level) - static_cast(min_level); - return *this; - } - - template - inline auto start(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 1 : 0; - } - int value = (((i - dec) >> m_shift) << m_shift) + dec; - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto end(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 0 : 1; - } - int value = (((i - dec) >> m_shift) + dec) << m_shift; - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto goback(const Func& f) const - { - auto new_f = [&, f](auto level, auto i) - { - auto [prev_lev, v] = f(level, i); - - int min_shift = static_cast(m_min_level) - static_cast(prev_lev); - int max_shift = static_cast(m_level) - static_cast(m_min_level); - - if constexpr (end) - { - i = end_shift(end_shift(v, min_shift), max_shift); - } - else - { - i = start_shift(start_shift(v, min_shift), max_shift); - } - - return std::make_pair(m_level, i); - }; - return new_f; - } - - std::size_t m_level; - int m_shift; - std::size_t m_min_level; - }; - - template - struct start_end_translate_function - { - using container_t = xt::xtensor_fixed>; - - explicit start_end_translate_function(const container_t& t) - : m_level(0) - , m_min_level(0) - , m_max_level(0) - , m_t(t) - { - } - - auto& operator()(auto level, auto min_level, auto max_level) - { - m_level = level; - m_min_level = min_level; - m_max_level = max_level; - return *this; - } - - template - inline auto start(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - int max2curr = static_cast(m_max_level) - static_cast(level); - int curr2min = static_cast(level) - static_cast(m_min_level); - int min2max = static_cast(m_max_level) - static_cast(m_min_level); - - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 1 : 0; - } - int value = (((((i - dec) >> max2curr) + m_t[d - 1]) >> curr2min) + dec) << min2max; - - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto end(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - int max2curr = static_cast(m_max_level) - static_cast(level); - int curr2min = static_cast(level) - static_cast(m_min_level); - int min2max = static_cast(m_max_level) - static_cast(m_min_level); - - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 0 : 1; - } - int value = (((((i - dec) >> max2curr) + m_t[d - 1]) >> curr2min) + dec) << min2max; - - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto goback(const Func& f) const - { - auto new_f = [&, f](auto level, auto i) - { - auto [prev_lev, v] = f(level, i); - - auto min_shift = static_cast(m_min_level) - static_cast(prev_lev); - auto max_shift = static_cast(m_level) - static_cast(m_min_level); - - if constexpr (end) - { - i = end_shift(end_shift(v, min_shift), max_shift) - m_t[d - 1]; - } - else - { - i = start_shift(start_shift(v, min_shift), max_shift) - m_t[d - 1]; - } - - return std::make_pair(m_level, i); - }; - return new_f; - } - - std::size_t m_level; - std::size_t m_min_level; - std::size_t m_max_level; - xt::xtensor_fixed> m_t; - }; - - template - struct start_end_contraction_function - { - explicit start_end_contraction_function(int c) - : m_level(0) - , m_min_level(0) - , m_max_level(0) - , m_c(c) - { - } - - auto& operator()(auto level, auto min_level, auto max_level) - { - m_level = level; - m_min_level = min_level; - m_max_level = max_level; - return *this; - } - - template - inline auto start(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - int max2curr = static_cast(m_max_level) - static_cast(level); - int curr2min = static_cast(level) - static_cast(m_min_level); - int min2max = static_cast(m_max_level) - static_cast(m_min_level); - - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 1 : 0; - dec = 1; - } - int value = (((((i - dec) >> max2curr) + m_c) >> curr2min) + dec) << min2max; - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto end(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - int max2curr = static_cast(m_max_level) - static_cast(level); - int curr2min = static_cast(level) - static_cast(m_min_level); - int min2max = static_cast(m_max_level) - static_cast(m_min_level); - - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 0 : 1; - } - int value = (((((i - dec) >> max2curr) - m_c) >> curr2min) + dec) << min2max; - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto goback(const Func& f) const - { - auto new_f = [&, f](auto level, auto i) - { - auto [prev_lev, v] = f(level, i); - - auto min_shift = static_cast(m_min_level) - static_cast(prev_lev); - auto max_shift = static_cast(m_level) - static_cast(m_min_level); - - if constexpr (end) - { - i = end_shift(end_shift(v, min_shift), max_shift) + m_c; - } - else - { - i = start_shift(start_shift(v, min_shift), max_shift) + m_c; - } - - return std::make_pair(m_level, i); - }; - return new_f; - } - - std::size_t m_level; - std::size_t m_min_level; - std::size_t m_max_level; - int m_c; - }; -} diff --git a/include/samurai/subset/translation.hpp b/include/samurai/subset/translation.hpp new file mode 100644 index 000000000..208774513 --- /dev/null +++ b/include/samurai/subset/translation.hpp @@ -0,0 +1,95 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_base.hpp" +#include "traversers/translation_traverser.hpp" + +#include +using namespace xt::placeholders; // this makes `_` available + +namespace samurai +{ + + template + class Translation; + + template + struct SetTraits> + { + static_assert(IsSet::value); + + template + using traverser_t = TranslationTraverser>; + + static constexpr std::size_t dim() + { + return Set::dim; + } + }; + + template + class Translation : public SetBase> + { + using Self = Translation; + + public: + + SAMURAI_SET_TYPEDEFS + + using translation_t = xt::xtensor_fixed>; + + template + Translation(const Set& set, const translation_expr_t& translation_expr) + : m_set(set) + , m_translation(translation_expr) + { + } + + inline std::size_t level_impl() const + { + return m_set.level(); + } + + inline bool exist_impl() const + { + return m_set.exist(); + } + + inline bool empty_impl() const + { + return m_set.empty(); + } + + template + inline void init_get_traverser_work_impl(const std::size_t n_traversers, std::integral_constant d_ic) const + { + m_set.init_get_traverser_work(n_traversers, d_ic); + } + + template + inline void clear_get_traverser_work_impl(std::integral_constant d_ic) const + { + m_set.clear_get_traverser_work(d_ic); + } + + template + inline traverser_t get_traverser_impl(const index_t& index, std::integral_constant d_ic) const + { + return traverser_t(m_set.get_traverser_impl(index - xt::view(m_translation, xt::range(1, _)), d_ic), m_translation[d]); + } + + private: + + Set m_set; + translation_t m_translation; + }; + + template + auto translate(const Set& set, const typename Translation>::translation_t& translation) + { + return Translation(self(set), translation); + } + +} // namespace samurai diff --git a/include/samurai/subset/traversers/box_traverser.hpp b/include/samurai/subset/traversers/box_traverser.hpp new file mode 100644 index 000000000..f74e2c4aa --- /dev/null +++ b/include/samurai/subset/traversers/box_traverser.hpp @@ -0,0 +1,63 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../../interval.hpp" +#include "set_traverser_base.hpp" +#include + +namespace samurai +{ + template + class Box; + + template + class BoxTraverser; + + template + struct SetTraverserTraits> + { + static_assert(std::same_as, B>); + + using interval_t = Interval; + using current_interval_t = const interval_t&; + }; + + template + class BoxTraverser : public SetTraverserBase> + { + using Self = BoxTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + BoxTraverser(const value_t& start, const value_t& end) + : m_current_interval{start, end} + , m_empty(start < end) + { + } + + inline bool is_empty_impl() const + { + return m_empty; + } + + inline void next_interval_impl() + { + m_empty = true; + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + interval_t m_current_interval; + bool m_empty; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/contraction_traverser.hpp b/include/samurai/subset/traversers/contraction_traverser.hpp new file mode 100644 index 000000000..8d9876f8d --- /dev/null +++ b/include/samurai/subset/traversers/contraction_traverser.hpp @@ -0,0 +1,71 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class ContractionTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverser::interval_t; + using current_interval_t = interval_t; + }; + + template + class ContractionTraverser : public SetTraverserBase> + { + using Self = ContractionTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + ContractionTraverser(const SetTraverser& set_traverser, const value_t contraction) + : m_set_traverser(set_traverser) + , m_contraction(contraction) + { + assert(m_contraction >= 0); + advance_to_next_valid_interval(); + } + + inline bool is_empty_impl() const + { + return m_set_traverser.is_empty(); + } + + inline void next_interval_impl() + { + m_set_traverser.next_interval(); + advance_to_next_valid_interval(); + } + + inline current_interval_t current_interval_impl() const + { + return current_interval_t(m_set_traverser.current_interval().start + m_contraction, + m_set_traverser.current_interval().end - m_contraction); + } + + private: + + inline void advance_to_next_valid_interval() + { + while (!m_set_traverser.is_empty() && m_set_traverser.current_interval().size() <= size_t(2 * m_contraction)) + { + m_set_traverser.next_interval(); + } + } + + SetTraverser m_set_traverser; + value_t m_contraction; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/difference_id_traverser.hpp b/include/samurai/subset/traversers/difference_id_traverser.hpp new file mode 100644 index 000000000..af2f9b294 --- /dev/null +++ b/include/samurai/subset/traversers/difference_id_traverser.hpp @@ -0,0 +1,66 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class DifferenceIdTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + static_assert((IsSetTraverser::value and ...)); + + using interval_t = typename FirstSetTraverser::interval_t; + using current_interval_t = interval_t; + }; + + template + class DifferenceIdTraverser : public SetTraverserBase> + { + using Self = DifferenceIdTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + static constexpr std::size_t nIntervals = 1 + sizeof...(OtherSetTraversers); + + DifferenceIdTraverser(const std::array& shifts, + const FirstSetTraverser& set_traverser, + const OtherSetTraversers&...) + : m_set_traverser(set_traverser) + , m_shift(shifts[0]) + { + } + + DifferenceIdTraverser() = delete; + + inline bool is_empty_impl() const + { + return m_set_traverser.is_empty(); + } + + inline void next_interval_impl() + { + m_set_traverser.next_interval(); + } + + inline current_interval_t current_interval_impl() const + { + return current_interval_t{m_set_traverser.current_interval().start << m_shift, m_set_traverser.current_interval().end << m_shift}; + } + + private: + + FirstSetTraverser m_set_traverser; + std::size_t m_shift; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/difference_traverser.hpp b/include/samurai/subset/traversers/difference_traverser.hpp new file mode 100644 index 000000000..6f3c7bdd3 --- /dev/null +++ b/include/samurai/subset/traversers/difference_traverser.hpp @@ -0,0 +1,128 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../../static_algorithm.hpp" +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class DifferenceTraverser; + + template + struct SetTraverserTraits> + { + static_assert((IsSetTraverser::value and ...)); + + using FirstSetTraverser = std::tuple_element_t<0, std::tuple>; + using interval_t = typename FirstSetTraverser::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class DifferenceTraverser : public SetTraverserBase> + { + using Self = DifferenceTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + using Childrens = std::tuple; + + template + using IthChild = std::tuple_element::type; + + static constexpr std::size_t nIntervals = std::tuple_size_v; + + DifferenceTraverser(const std::array& shifts, const SetTraversers&... set_traversers) + : m_min_start(std::numeric_limits::min()) + , m_set_traversers(set_traversers...) + , m_shifts(shifts) + { + compute_current_interval(); + } + + inline bool is_empty_impl() const + { + return std::get<0>(m_set_traversers).is_empty(); + } + + inline void next_interval_impl() + { + advance_ref_interval(); + compute_current_interval(); + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + inline void advance_ref_interval() + { + if (m_current_interval.end != std::get<0>(m_set_traversers).current_interval().end << m_shifts[0]) + { + // we have removed the beginning of the current interval. + // so ve remove [m_current_interval.start, m_current_interval.end) from std::get<0>(m_set_traversers).current_interval() + m_min_start = m_current_interval.end; + } + else + { + // all of the current interval has been removed. + // move to the next one. + std::get<0>(m_set_traversers).next_interval(); + } + } + + inline void compute_current_interval() + { + while (!std::get<0>(m_set_traversers).is_empty() && !try_to_compute_current_interval()) + { + advance_ref_interval(); + } + } + + inline bool try_to_compute_current_interval() + { + assert(!std::get<0>(m_set_traversers).is_empty()); + + m_current_interval.start = std::max(m_min_start, std::get<0>(m_set_traversers).current_interval().start << m_shifts[0]); + m_current_interval.end = std::get<0>(m_set_traversers).current_interval().end << m_shifts[0]; + + static_for<1, nIntervals>::apply( + [this](const auto i) + { + IthChild& set_traverser = std::get(m_set_traversers); + + while (!set_traverser.is_empty() && (set_traverser.current_interval().end << m_shifts[i]) < m_current_interval.start) + { + set_traverser.next_interval(); + } + + if (!set_traverser.is_empty() && (set_traverser.current_interval().start << m_shifts[i]) <= m_current_interval.start) + { + m_current_interval.start = set_traverser.current_interval().end << m_shifts[i]; + m_min_start = m_current_interval.start; + set_traverser.next_interval(); + } + if (!set_traverser.is_empty() && (set_traverser.current_interval().start << m_shifts[i]) <= m_current_interval.end) + { + m_current_interval.end = set_traverser.current_interval().start << m_shifts[i]; + } + }); + + return m_current_interval.is_valid(); + } + + interval_t m_current_interval; + value_t m_min_start; + Childrens m_set_traversers; + const std::array& m_shifts; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/expansion_traverser.hpp b/include/samurai/subset/traversers/expansion_traverser.hpp new file mode 100644 index 000000000..947369935 --- /dev/null +++ b/include/samurai/subset/traversers/expansion_traverser.hpp @@ -0,0 +1,100 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class ExpansionTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverserTraits::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class ExpansionTraverser : public SetTraverserBase> + { + using Self = ExpansionTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + using SetTraverserIterator = typename std::vector::iterator; + + ExpansionTraverser(SetTraverserIterator begin_set_traverser, SetTraverserIterator end_set_traverser, const value_t expansion) + : m_set_traversers(begin_set_traverser, end_set_traverser) + , m_expansion(expansion) + { + assert(m_expansion > 0); + next_interval_impl(); + } + + inline bool is_empty_impl() const + { + return m_current_interval.start == std::numeric_limits::max(); + } + + inline void next_interval_impl() + { + m_current_interval.start = std::numeric_limits::max(); + + // We find the start of the interval, i.e. the smallest set_traverser.current_interval().start + for (const SetTraverser& set_traverser : m_set_traversers) + { + if (!set_traverser.is_empty() && (set_traverser.current_interval().start - m_expansion < m_current_interval.start)) + { + m_current_interval.start = set_traverser.current_interval().start - m_expansion; + m_current_interval.end = set_traverser.current_interval().end + m_expansion; + } + } + // Now we find the end of the interval, i.e. the largest set_traverser.current_interval().end + // such that set_traverser.current_interval().start - expansion < m_current_interval.end + bool is_done = false; + while (!is_done) + { + is_done = true; + // advance set traverses that are behind current interval + for (SetTraverser& set_traverser : m_set_traversers) + { + while (!set_traverser.is_empty() && set_traverser.current_interval().end + m_expansion <= m_current_interval.end) + { + set_traverser.next_interval(); + } + } + // try to find a new end + for (const SetTraverser& set_traverser : m_set_traversers) + { + // there is an overlap + if (!set_traverser.is_empty() && set_traverser.current_interval().start - m_expansion <= m_current_interval.end) + { + is_done = false; + m_current_interval.end = set_traverser.current_interval().end + m_expansion; + } + } + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + std::span m_set_traversers; + value_t m_expansion; + interval_t m_current_interval; + bool m_isEmpty; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/intersection_traverser.hpp b/include/samurai/subset/traversers/intersection_traverser.hpp new file mode 100644 index 000000000..54c9a3089 --- /dev/null +++ b/include/samurai/subset/traversers/intersection_traverser.hpp @@ -0,0 +1,107 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class IntersectionTraverser; + + template + struct SetTraverserTraits> + { + static_assert((IsSetTraverser::value and ...)); + + using FirstSetTraverser = std::tuple_element_t<0, std::tuple>; + using interval_t = typename FirstSetTraverser::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class IntersectionTraverser : public SetTraverserBase> + { + using Self = IntersectionTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + using Childrens = std::tuple; + + template + using IthChild = std::tuple_element::type; + + static constexpr std::size_t nIntervals = std::tuple_size_v; + + IntersectionTraverser(const std::array& shifts, const SetTraversers&... set_traverser) + : m_set_traversers(set_traverser...) + , m_shifts(shifts) + { + next_interval_impl(); + } + + inline bool is_empty_impl() const + { + return !m_current_interval.is_valid(); + } + + inline void next_interval_impl() + { + m_current_interval.start = 0; + m_current_interval.end = 0; + + while (not_is_any_child_empty() && m_current_interval.start >= m_current_interval.end) + { + m_current_interval.start = std::numeric_limits::min(); + m_current_interval.end = std::numeric_limits::max(); + + enumerate_const_items( + m_set_traversers, + [this](const auto i, const auto& set_traverser) + { + if (!set_traverser.is_empty()) + { + m_current_interval.start = std::max(m_current_interval.start, + set_traverser.current_interval().start << m_shifts[i]); + m_current_interval.end = std::min(m_current_interval.end, set_traverser.current_interval().end << m_shifts[i]); + } + }); + + enumerate_items( + m_set_traversers, + [this](const auto i, auto& set_traverser) + { + if (!set_traverser.is_empty() && ((set_traverser.current_interval().end << m_shifts[i]) == m_current_interval.end)) + { + set_traverser.next_interval(); + } + }); + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + inline bool not_is_any_child_empty() const + { + return std::apply( + [](const auto&... set_traversers) + { + return (!set_traversers.is_empty() and ...); + }, + m_set_traversers); + } + + interval_t m_current_interval; + Childrens m_set_traversers; + const std::array& m_shifts; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/last_dim_expansion_traverser.hpp b/include/samurai/subset/traversers/last_dim_expansion_traverser.hpp new file mode 100644 index 000000000..ab47787a7 --- /dev/null +++ b/include/samurai/subset/traversers/last_dim_expansion_traverser.hpp @@ -0,0 +1,76 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class LastDimExpansionTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverserTraits::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class LastDimExpansionTraverser : public SetTraverserBase> + { + using Self = LastDimExpansionTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + LastDimExpansionTraverser(const SetTraverser& set_traverser, const value_t expansion) + : m_set_traverser(set_traverser) + , m_expansion(expansion) + { + assert(m_expansion > 0); + next_interval_impl(); + } + + inline bool is_empty_impl() const + { + return m_isEmpty; + } + + inline void next_interval_impl() + { + m_isEmpty = m_set_traverser.is_empty(); + + if (!m_isEmpty) + { + m_current_interval.start = m_set_traverser.current_interval().start - m_expansion; + m_current_interval.end = m_set_traverser.current_interval().end + m_expansion; + + m_set_traverser.next_interval(); + while (!m_set_traverser.is_empty() and m_set_traverser.current_interval().start - m_expansion <= m_current_interval.end) + { + m_current_interval.end = m_set_traverser.current_interval().end + m_expansion; + m_set_traverser.next_interval(); + } + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + SetTraverser m_set_traverser; + value_t m_expansion; + interval_t m_current_interval; + bool m_isEmpty; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/last_dim_projection_loi_traverser.hpp b/include/samurai/subset/traversers/last_dim_projection_loi_traverser.hpp new file mode 100644 index 000000000..07847cc03 --- /dev/null +++ b/include/samurai/subset/traversers/last_dim_projection_loi_traverser.hpp @@ -0,0 +1,118 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + template + class LastDimProjectionLOITraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverser::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class LastDimProjectionLOITraverser : public SetTraverserBase> + { + using Self = LastDimProjectionLOITraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + LastDimProjectionLOITraverser(const SetTraverser& set_traverser, const ProjectionType projectionType, const std::size_t shift) + : m_set_traverser(set_traverser) + , m_projectionType(projectionType) + , m_shift(shift) + , m_isEmpty(set_traverser.is_empty()) + { + if (m_projectionType == ProjectionType::COARSEN) + { + next_interval_coarsen(); + } + else if (!m_isEmpty) + { + m_current_interval.start = m_set_traverser.current_interval().start << shift; + m_current_interval.end = m_set_traverser.current_interval().end << shift; + } + } + + inline bool is_empty_impl() const + { + return m_isEmpty; + } + + inline void next_interval_impl() + { + if (m_projectionType == ProjectionType::COARSEN) + { + next_interval_coarsen(); + } + else + { + m_set_traverser.next_interval(); + m_isEmpty = m_set_traverser.is_empty(); + if (!m_isEmpty) + { + m_current_interval.start = m_set_traverser.current_interval().start << m_shift; + m_current_interval.end = m_set_traverser.current_interval().end << m_shift; + } + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + inline void next_interval_coarsen() + { + if (!m_set_traverser.is_empty()) + { + m_current_interval.start = coarsen_start(m_set_traverser.current_interval()); + m_current_interval.end = coarsen_end(m_set_traverser.current_interval()); + + m_set_traverser.next_interval(); + + // when coarsening, two disjoint intervals may be merged. + // we need to check if the next_interval overlaps + for (; !m_set_traverser.is_empty() && coarsen_start(m_set_traverser.current_interval()) <= m_current_interval.end; + m_set_traverser.next_interval()) + { + m_current_interval.end = coarsen_end(m_set_traverser.current_interval()); + } + } + else + { + m_isEmpty = true; + } + } + + inline value_t coarsen_start(const interval_t& interval) const + { + return interval.start >> m_shift; + } + + inline value_t coarsen_end(const interval_t& interval) const + { + return ((interval.end - 1) >> m_shift) + 1; + } + + SetTraverser m_set_traverser; + ProjectionType m_projectionType; + std::size_t m_shift; + interval_t m_current_interval; + bool m_isEmpty; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/lca_traverser.hpp b/include/samurai/subset/traversers/lca_traverser.hpp new file mode 100644 index 000000000..00a1f3620 --- /dev/null +++ b/include/samurai/subset/traversers/lca_traverser.hpp @@ -0,0 +1,62 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + template + class LevelCellArray; + + template + class LCATraverser; + + template + struct SetTraverserTraits> + { + static_assert(std::same_as, LCA>); + + using interval_t = typename LCA::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class LCATraverser : public SetTraverserBase> + { + using Self = LCATraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + using interval_iterator = typename std::vector::const_iterator; + + LCATraverser(const interval_iterator first, const interval_iterator end) + : m_first_interval(first) + , m_end_interval(end) + { + } + + inline bool is_empty_impl() const + { + return m_first_interval == m_end_interval; + } + + inline void next_interval_impl() + { + ++m_first_interval; + } + + inline current_interval_t current_interval_impl() const + { + return *m_first_interval; + } + + private: + + interval_iterator m_first_interval; + interval_iterator m_end_interval; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/projection_loi_traverser.hpp b/include/samurai/subset/traversers/projection_loi_traverser.hpp new file mode 100644 index 000000000..32afbc79f --- /dev/null +++ b/include/samurai/subset/traversers/projection_loi_traverser.hpp @@ -0,0 +1,87 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../../list_of_intervals.hpp" +#include "set_traverser_base.hpp" + +namespace samurai +{ + template + class ProjectionLOITraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverser::interval_t; + using current_interval_t = interval_t; + }; + + template + class ProjectionLOITraverser : public SetTraverserBase> + { + using Self = ProjectionLOITraverser; + using SetTraverserIterator = typename ListOfIntervals::const_iterator; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + ProjectionLOITraverser(SetTraverser set_traverser, const ProjectionType projectionType, const std::size_t shift) + : m_set_traverser(set_traverser) + , m_shift(shift) + , m_projectionType(ProjectionType::REFINE) + { + assert(projectionType == m_projectionType); + } + + ProjectionLOITraverser(SetTraverser set_traverser, SetTraverserIterator first_interval, SetTraverserIterator bound_interval) + : m_set_traverser(set_traverser) + , m_first_interval(first_interval) + , m_bound_interval(bound_interval) + , m_projectionType(ProjectionType::COARSEN) + { + } + + inline bool is_empty_impl() const + { + if (m_projectionType == ProjectionType::COARSEN) + { + return m_first_interval == m_bound_interval; + } + else + { + return m_set_traverser.is_empty(); + } + } + + inline void next_interval_impl() + { + if (m_projectionType == ProjectionType::COARSEN) + { + ++m_first_interval; + } + else + { + m_set_traverser.next_interval(); + } + } + + inline current_interval_t current_interval_impl() const + { + return (m_projectionType == ProjectionType::COARSEN) ? *m_first_interval : m_set_traverser.current_interval() << m_shift; + } + + private: + + SetTraverser m_set_traverser; // only used when refining + std::size_t m_shift; // only used when refining + SetTraverserIterator m_first_interval; // only use when coarsening + SetTraverserIterator m_bound_interval; // only use when coarsening + ProjectionType m_projectionType; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/projection_traverser.hpp b/include/samurai/subset/traversers/projection_traverser.hpp new file mode 100644 index 000000000..267e0df1e --- /dev/null +++ b/include/samurai/subset/traversers/projection_traverser.hpp @@ -0,0 +1,162 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../projection_type.hpp" +#include "set_traverser_base.hpp" + +namespace samurai +{ + template + class ProjectionTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverser::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class ProjectionTraverser : public SetTraverserBase> + { + using Self = ProjectionTraverser; + using SetTraverserIterator = typename std::vector::iterator; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + ProjectionTraverser(SetTraverserIterator set_traverser, const ProjectionType projectionType, const std::size_t shift) + : m_set_traversers(set_traverser, set_traverser + 1) + , m_projectionType(projectionType) + , m_shift(shift) + , m_isEmpty(set_traverser->is_empty()) + { + if (!m_isEmpty) + { + if (m_projectionType == ProjectionType::COARSEN) + { + m_current_interval.start = coarsen_start(m_set_traversers[0].current_interval()); + m_current_interval.end = coarsen_end(m_set_traversers[0].current_interval()); + + m_set_traversers[0].next_interval(); + + // when coarsening, two disjoint intervals may be merged. + // we need to check if the next_interval overlaps + for (; !m_set_traversers[0].is_empty() && coarsen_start(m_set_traversers[0].current_interval()) <= m_current_interval.end; + m_set_traversers[0].next_interval()) + { + m_current_interval.end = coarsen_end(m_set_traversers[0].current_interval()); + } + } + else + { + m_current_interval.start = m_set_traversers[0].current_interval().start << shift; + m_current_interval.end = m_set_traversers[0].current_interval().end << shift; + } + } + } + + /* + * This constructor only works for coarsening + */ + ProjectionTraverser(SetTraverserIterator begin_set_traversers, SetTraverserIterator end_set_traversers, const std::size_t shift) + : m_set_traversers(begin_set_traversers, end_set_traversers) + , m_projectionType(ProjectionType::COARSEN) + , m_shift(shift) + { + next_interval_coarsen(); + } + + inline bool is_empty_impl() const + { + return m_isEmpty; + } + + inline void next_interval_impl() + { + if (m_projectionType == ProjectionType::COARSEN) + { + next_interval_coarsen(); + } + else + { + m_set_traversers[0].next_interval(); + m_isEmpty = m_set_traversers[0].is_empty(); + if (!m_isEmpty) + { + m_current_interval.start = m_set_traversers[0].current_interval().start << m_shift; + m_current_interval.end = m_set_traversers[0].current_interval().end << m_shift; + } + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + inline void next_interval_coarsen() + { + m_current_interval.start = std::numeric_limits::max(); + // We find the start of the interval, i.e. the smallest set_traverser.current_interval().start >> m_shift + for (const SetTraverser& set_traverser : m_set_traversers) + { + if (!set_traverser.is_empty() && (coarsen_start(set_traverser.current_interval()) < m_current_interval.start)) + { + m_current_interval.start = coarsen_start(set_traverser.current_interval()); + m_current_interval.end = coarsen_end(set_traverser.current_interval()); + } + } + // Now we find the end of the interval, i.e. the largest set_traverser.current_interval().end >> m_shift + // such that (set_traverser.current_interval().start >> m_shift) < m_current_interval.end + bool is_done = false; + while (!is_done) + { + is_done = true; + // advance set traverses that are behind current interval + for (SetTraverser& set_traverser : m_set_traversers) + { + while (!set_traverser.is_empty() && (coarsen_end(set_traverser.current_interval()) <= m_current_interval.end)) + { + set_traverser.next_interval(); + } + } + // try to find a new end + for (const SetTraverser& set_traverser : m_set_traversers) + { + // there is an overlap + if (!set_traverser.is_empty() && (coarsen_start(set_traverser.current_interval()) <= m_current_interval.end)) + { + is_done = false; + m_current_interval.end = coarsen_end(set_traverser.current_interval()); + } + } + } + m_isEmpty = (m_current_interval.start == std::numeric_limits::max()); + } + + inline value_t coarsen_start(const interval_t& interval) const + { + return interval.start >> m_shift; + } + + inline value_t coarsen_end(const interval_t& interval) const + { + return ((interval.end - 1) >> m_shift) + 1; + } + + std::span m_set_traversers; + ProjectionType m_projectionType; + std::size_t m_shift; + interval_t m_current_interval; + bool m_isEmpty; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/set_traverser_base.hpp b/include/samurai/subset/traversers/set_traverser_base.hpp new file mode 100644 index 000000000..2d7855293 --- /dev/null +++ b/include/samurai/subset/traversers/set_traverser_base.hpp @@ -0,0 +1,61 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include +#include + +namespace samurai +{ + template + struct SetTraverserTraits; + + template + class SetTraverserBase + { + public: + + using interval_t = typename SetTraverserTraits::interval_t; + using current_interval_t = typename SetTraverserTraits::current_interval_t; + using value_t = typename interval_t::value_t; + + const Derived& derived_cast() const + { + return static_cast(*this); + } + + Derived& derived_cast() + { + return static_cast(*this); + } + + inline bool is_empty() const + { + return derived_cast().is_empty_impl(); + } + + inline void next_interval() + { + assert(!is_empty()); + derived_cast().next_interval_impl(); + } + + inline current_interval_t current_interval() const + { + return derived_cast().current_interval_impl(); + } + }; + + template + struct IsSetTraverser : std::bool_constant, T>::value> + { + }; + +#define SAMURAI_SET_TRAVERSER_TYPEDEFS \ + using Base = SetTraverserBase; \ + using interval_t = typename Base::interval_t; \ + using current_interval_t = typename Base::current_interval_t; \ + using value_t = typename Base::value_t; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/translation_traverser.hpp b/include/samurai/subset/traversers/translation_traverser.hpp new file mode 100644 index 000000000..3f7b6cbb4 --- /dev/null +++ b/include/samurai/subset/traversers/translation_traverser.hpp @@ -0,0 +1,60 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class TranslationTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverserTraits::interval_t; + using current_interval_t = interval_t; + }; + + template + class TranslationTraverser : public SetTraverserBase> + { + using Self = TranslationTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + TranslationTraverser(const SetTraverser& set_traverser, const value_t& translation) + : m_set_traverser(set_traverser) + , m_translation(translation) + { + } + + inline bool is_empty_impl() const + { + return m_set_traverser.is_empty(); + } + + inline void next_interval_impl() + { + m_set_traverser.next_interval(); + } + + inline current_interval_t current_interval_impl() const + { + return current_interval_t{m_set_traverser.current_interval().start + m_translation, + m_set_traverser.current_interval().end + m_translation}; + } + + private: + + SetTraverser m_set_traverser; + value_t m_translation; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/union_traverser.hpp b/include/samurai/subset/traversers/union_traverser.hpp new file mode 100644 index 000000000..00dc22b81 --- /dev/null +++ b/include/samurai/subset/traversers/union_traverser.hpp @@ -0,0 +1,108 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class UnionTraverser; + + template + struct SetTraverserTraits> + { + static_assert((IsSetTraverser::value and ...)); + + using FirstSetTraverser = std::tuple_element_t<0, std::tuple>; + using interval_t = typename FirstSetTraverser::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class UnionTraverser : public SetTraverserBase> + { + using Self = UnionTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + using Childrens = std::tuple; + + template + using IthChild = std::tuple_element::type; + + static constexpr std::size_t nIntervals = std::tuple_size::value; + + UnionTraverser(const std::array& shifts, const SetTraversers&... set_traversers) + : m_set_traversers(set_traversers...) + , m_shifts(shifts) + { + next_interval_impl(); + } + + inline bool is_empty_impl() const + { + return m_current_interval.start == std::numeric_limits::max(); + } + + inline void next_interval_impl() + { + m_current_interval.start = std::numeric_limits::max(); + // We find the start of the interval, i.e. the smallest set_traverser.current_interval().start << m_shifts[i] + enumerate_const_items( + m_set_traversers, + [this](const auto i, const auto& set_traverser) + { + if (!set_traverser.is_empty() && ((set_traverser.current_interval().start << m_shifts[i]) < m_current_interval.start)) + { + m_current_interval.start = set_traverser.current_interval().start << m_shifts[i]; + m_current_interval.end = set_traverser.current_interval().end << m_shifts[i]; + } + }); + // Now we find the end of the interval, i.e. the largest set_traverser.current_interval().end << m_shifts[i] + // such that (set_traverser.current_interval().start << m_shifts[i]) < m_current_interval.end + bool is_done = false; + while (!is_done) + { + is_done = true; + // advance set traverses that are behind current interval + enumerate_items( + m_set_traversers, + [this](const auto i, auto& set_traverser) + { + while (!set_traverser.is_empty() && (set_traverser.current_interval().end << m_shifts[i]) <= m_current_interval.end) + { + set_traverser.next_interval(); + } + }); + // try to find a new end + enumerate_const_items( + m_set_traversers, + [&is_done, this](const auto i, const auto& set_traverser) + { + // there is an overlap + if (!set_traverser.is_empty() && (set_traverser.current_interval().start << m_shifts[i]) <= m_current_interval.end) + { + is_done = false; + m_current_interval.end = set_traverser.current_interval().end << m_shifts[i]; + } + }); + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + interval_t m_current_interval; + Childrens m_set_traversers; + const std::array& m_shifts; + }; + +} // namespace samurai diff --git a/include/samurai/subset/utils.hpp b/include/samurai/subset/utils.hpp index 4c774c8ef..a45c1e1bb 100644 --- a/include/samurai/subset/utils.hpp +++ b/include/samurai/subset/utils.hpp @@ -3,73 +3,55 @@ #pragma once -#include -#include -#include -#include +#include +#include -// namespace samurai::experimental namespace samurai { - template - static constexpr T sentinel = std::numeric_limits::max(); - - template - inline T end_shift(T value, T shift) + //////////////////////////////////////////////////////////////////////// + //// misc + //////////////////////////////////////////////////////////////////////// + template + const T& vmin(const T& a) { - return shift >= 0 ? value << shift : ((value - 1) >> -shift) + 1; + return a; } - template - inline T start_shift(T value, T shift) + template + const T& vmax(const T& a) { - return shift >= 0 ? value << shift : value >> -shift; + return a; } - constexpr auto compute_min(auto const& value, auto const&... args) + template + T vmin(const T& a, const T& b, const Ts&... others) { - if constexpr (sizeof...(args) == 0u) // Single argument case! + if constexpr (sizeof...(Ts) == 0u) { - return value; + return std::min(a, b); } - else // For the Ts... + else { - const auto min = compute_min(args...); - return value < min ? value : min; + return a < b ? vmin(a, others...) : vmin(b, others...); } } - constexpr auto compute_max(auto const& value, auto const&... args) + template + T vmax(const T& a, const T& b, const Ts&... others) { - if constexpr (sizeof...(args) == 0u) // Single argument case! + if constexpr (sizeof...(Ts) == 0u) { - return value; + return std::max(a, b); } - else // For the Ts... + else { - const auto max = compute_max(args...); - return value > max ? value : max; + return a > b ? vmax(a, others...) : vmax(b, others...); } } - template - struct get_interval_type - { - using type = typename S1::interval_t; - }; - - template - using get_interval_t = typename get_interval_type::type; - - template - struct get_set_dim - { - static constexpr std::size_t value = S1::dim; - }; - - template - constexpr std::size_t get_set_dim_v = get_set_dim...>::value; - + //////////////////////////////////////////////////////////////////////// + //// intervals args + //////////////////////////////////////////////////////////////////////// template ::value_type::value_type> ForwardIt lower_bound_interval(ForwardIt begin, ForwardIt end, const T& value) { @@ -94,20 +76,40 @@ namespace samurai }); } - template - constexpr void zip_apply_impl(F&& f, Tuple1&& t1, Tuple2&& t2, std::index_sequence) + //////////////////////////////////////////////////////////////////////// + //// tuple iteration + //////////////////////////////////////////////////////////////////////// + namespace detail { - (f(std::get(std::forward(t1)), std::get(std::forward(t2))), ...); + template + Func enumerate_items(Tuple& tuple, Func func, std::index_sequence) + { + (func(Is, std::get(tuple)), ...); + return func; + } + + template + Func enumerate_const_items(const Tuple& tuple, Func func, std::index_sequence) + { + (func(Is, std::get(tuple)), ...); + return func; + } } - template - constexpr void zip_apply(F&& f, Tuple1&& t1, Tuple2&& t2) + template + Func enumerate_items(Tuple& tuple, Func func) { - constexpr std::size_t size1 = std::tuple_size_v>; - constexpr std::size_t size2 = std::tuple_size_v>; + constexpr std::size_t N = std::tuple_size_v>; - static_assert(size1 == size2, "Tuples must have the same size"); + return detail::enumerate_items(tuple, func, std::make_index_sequence{}); + } + + template + Func enumerate_const_items(const Tuple& tuple, Func func) + { + constexpr std::size_t N = std::tuple_size_v>; - zip_apply_impl(std::forward(f), std::forward(t1), std::forward(t2), std::make_index_sequence{}); + return detail::enumerate_const_items(tuple, func, std::make_index_sequence{}); } -} + +} // namespace samurai diff --git a/include/samurai/subset/visitor.hpp b/include/samurai/subset/visitor.hpp deleted file mode 100644 index 507c59f3f..000000000 --- a/include/samurai/subset/visitor.hpp +++ /dev/null @@ -1,379 +0,0 @@ -// Copyright 2018-2025 the samurai's authors -// SPDX-License-Identifier: BSD-3-Clause - -#pragma once - -#include -#include - -#include "utils.hpp" - -namespace samurai -{ - template - class IntervalListRange - { - public: - - using container_t = container_; - using value_t = typename container_t::value_type; - using iterator_t = typename container_t::const_iterator; - - IntervalListRange(const container_t& data, std::ptrdiff_t start, std::ptrdiff_t end) - : m_data(data) - , m_work(data) - , m_begin(data.cbegin() + start) - , m_end(data.cbegin() + end) - { - } - - IntervalListRange(const container_t& data, const container_t& w) - : m_data(data) - , m_work(w) - , m_begin(w.cbegin()) - , m_end(w.cend()) - { - } - - inline auto begin() const - { - return m_begin; - } - - inline auto end() const - { - return m_end; - } - - private: - - const container_t& m_data; - const container_t& m_work; - iterator_t m_begin; - iterator_t m_end; - }; - - template - class IntervalListVisitor - { - public: - - using container_t = container_; - using base_t = IntervalListRange; - using iterator_t = typename base_t::iterator_t; - using interval_t = typename base_t::value_t; - using value_t = typename interval_t::value_t; - - IntervalListVisitor(auto lca_level, auto level, auto max_level, const IntervalListRange& intervals) - : m_lca_level(static_cast(lca_level)) - , m_shift2dest(static_cast(max_level) - static_cast(level)) - , m_shift2ref(static_cast(max_level) - static_cast(lca_level)) - , m_intervals(intervals) - , m_first(intervals.begin()) - , m_last(intervals.end()) - , m_current(std::numeric_limits::min()) - , m_is_start(true) - { - } - - explicit IntervalListVisitor(IntervalListRange&& intervals) - : m_lca_level(std::numeric_limits::infinity()) - , m_shift2dest(std::numeric_limits::infinity()) - , m_shift2ref(std::numeric_limits::infinity()) - , m_intervals(std::move(intervals)) - , m_first(m_intervals.begin()) - , m_last(m_intervals.end()) - , m_current(sentinel) - , m_is_start(true) - { - } - - template - inline auto start(const auto& it, Func& start_fct) const - { - auto i = it->start << m_shift2ref; - return start_fct(m_lca_level, i, 0); - } - - template - inline auto end(const auto& it, Func& end_fct) const - { - auto i = it->end << m_shift2ref; - return end_fct(m_lca_level, i, 1); - } - - inline bool is_in(auto scan) const - { - // Recall that we check if scan is inside an interval defined as [start, - // end[. The end of the interval is not included. - // - // if the m_current value is the start of the interval which means m_is_start = - // true then if scan is lower than m_current, scan is not in the - // interval. - // - // if the m_current value is the end of the interval which means m_is_start = false - // then if scan is lower than m_current, scan is in the interval. - return m_current != sentinel && !((scan < m_current) ^ (!m_is_start)); - } - - inline bool is_empty() const - { - return m_current == sentinel; - } - - inline auto min() const - { - return m_current; - } - - inline auto shift() const - { - return m_shift2dest; - } - - template - inline void next_interval(StartEnd& start_and_stop) - { - auto& [start_fct, end_fct] = start_and_stop; // cppcheck-suppress variableScope - - auto i_start = start(m_first, start_fct); - auto i_end = end(m_first, end_fct); - while (m_first + 1 != m_last && i_end >= start(m_first + 1, start_fct)) - { - ++m_first; - i_end = end(m_first, end_fct); - } - m_current_interval = {i_start, i_end}; - - if (m_current_interval.is_valid()) - { - m_current = m_current_interval.start; - } - else - { - m_current = sentinel; - } - } - - template - inline void next(auto scan, StartEnd& start_and_stop) - { - if (m_current == std::numeric_limits::min()) - { - next_interval(start_and_stop); - return; - } - - if (m_current == scan) - { - if (m_is_start) - { - m_current = m_current_interval.end; - } - else - { - ++m_first; - - if (m_first == m_last) - { - m_current = sentinel; - return; - } - next_interval(start_and_stop); - } - m_is_start = !m_is_start; - } - } - - private: - - int m_lca_level; - int m_shift2dest; - int m_shift2ref; - IntervalListRange m_intervals; - iterator_t m_first; - iterator_t m_last; - value_t m_current; - interval_t m_current_interval; - bool m_is_start; - }; - - template - class SetTraverser - { - public: - - static constexpr std::size_t dim = get_set_dim_v; - using set_type = std::tuple; - using interval_t = get_interval_t; - - SetTraverser(int shift, const Operator& op, S&&... s) - : m_shift(shift) - , m_operator(op) - , m_s(std::forward(s)...) - { - } - - inline auto shift() const - { - return m_shift; - } - - inline bool is_in(auto scan) const - { - return std::apply( - [this, scan](auto&&... args) - { - return m_operator.is_in(scan, args...); - }, - m_s); - } - - inline bool is_empty() const - { - return std::apply( - [this](auto&&... args) - { - return m_operator.is_empty(args...); - }, - m_s); - } - - inline auto min() const - { - return std::apply( - [](auto&&... args) - { - return compute_min(args.min()...); - }, - m_s); - } - - template - void next(auto scan, StartEnd&& start_and_stop) - { - zip_apply( - [scan](auto& arg, auto& start_end_fct) - { - arg.next(scan, start_end_fct); - }, - m_s, - std::forward(start_and_stop)); - } - - private: - - int m_shift; - Operator m_operator; - set_type m_s; - }; - - struct IntersectionOp - { - bool is_in(auto scan, const auto&... args) const - { - return (args.is_in(scan) && ...); - } - - bool is_empty(const auto&... args) const - { - return (args.is_empty() || ...); - } - - bool exist(const auto&... args) const - { - return (args.exist() && ...); - } - }; - - struct UnionOp - { - bool is_in(auto scan, const auto&... args) const - { - return (args.is_in(scan) || ...); - } - - bool is_empty(const auto&... args) const - { - return (args.is_empty() && ...); - } - - bool exist(const auto&... args) const - { - return (args.exist() || ...); - } - }; - - struct DifferenceOp - { - bool is_in(auto scan, const auto& arg, const auto&... args) const - { - return arg.is_in(scan) && !(args.is_in(scan) || ...); - } - - bool is_empty(const auto& arg, const auto&...) const - { - return arg.is_empty(); - } - - bool exist(const auto& arg, const auto&...) const - { - return arg.exist(); - } - }; - - struct Difference2Op - { - bool is_in(auto scan, const auto& arg, const auto&...) const - { - return arg.is_in(scan); - } - - bool is_empty(const auto& arg, const auto&...) const - { - return arg.is_empty(); - } - - bool exist(const auto& arg, const auto&...) const - { - return arg.exist(); - } - }; - - template - auto get_operator(const operator_t& op) - { - return op; - } - - template - auto get_operator(const DifferenceOp& op) - { - if constexpr (d == 1) - { - return op; - } - else - { - return Difference2Op(); - } - } - - struct SelfOp - { - bool is_in(auto scan, const auto& arg) const - { - return arg.is_in(scan); - } - - bool is_empty(const auto& arg) const - { - return arg.is_empty(); - } - - bool exist(const auto& arg) const - { - return arg.exist(); - } - }; -} diff --git a/tests/test_subset.cpp b/tests/test_subset.cpp index 1528938d4..e1f0c256a 100644 --- a/tests/test_subset.cpp +++ b/tests/test_subset.cpp @@ -14,6 +14,9 @@ #include #include +#include +#include + namespace samurai { TEST(subset, lower_bound) @@ -144,9 +147,9 @@ namespace samurai TEST(subset, compute_min) { - EXPECT_EQ(1, compute_min(3, 4, 1, 4)); - EXPECT_EQ(0, compute_min(0, 0, 0, 0)); - EXPECT_EQ(-1, compute_min(-1, -1, -1, -1)); + EXPECT_EQ(1, vmin(3, 4, 1, 4)); + EXPECT_EQ(0, vmin(0, 0, 0, 0)); + EXPECT_EQ(-1, vmin(-1, -1, -1, -1)); } TEST(subset, check_dim) @@ -224,8 +227,9 @@ namespace samurai EXPECT_EQ(interval_t(0, 20), i); }); - EXPECT_EQ(set.on(5).level(), 5); - apply(set, + auto set2 = set.on(5); + EXPECT_EQ(set2.level(), 5); + apply(set2, [](auto& i, auto) { EXPECT_EQ(interval_t(0, 40), i); @@ -454,6 +458,7 @@ namespace samurai {5, {5, 33}} }; std::size_t ie = 0; + fmt::print("====================================================\n"); apply(union_(ca_1[7], ca_2[7]).on(6), [&](auto& i, auto& index) { @@ -511,6 +516,134 @@ namespace samurai } } + TEST(subset, expand) + { + using interval_t = typename LevelCellArray<2>::interval_t; + using expected_t = std::vector>; + + LevelCellArray<2> ca; + + ca.add_interval_back({0, 1}, {0}); + + { + const auto translated_ca = translate(ca, {3 + 1, 0}); + const auto joined_cas = union_(ca, translated_ca); + + const auto set = expand(joined_cas, 3); + + expected_t expected{ + {-3, {-3, 8}}, + {-2, {-3, 8}}, + {-1, {-3, 8}}, + {0, {-3, 8}}, + {1, {-3, 8}}, + {2, {-3, 8}}, + {3, {-3, 8}} + }; + + bool is_set_empty = true; + std::size_t ie = 0; + set( + [&expected, &is_set_empty, &ie](const auto& x_interval, const auto& yz) + { + is_set_empty = false; + EXPECT_EQ(expected[ie++], std::make_pair(yz[0], x_interval)); + }); + EXPECT_EQ(ie, expected.size()); + EXPECT_FALSE(is_set_empty); + } + + { + const auto translated_ca = translate(ca, {0, 3 + 1}); + const auto joined_cas = union_(ca, translated_ca); + + const auto set = expand(joined_cas, 3); + + expected_t expected{ + {-3, {-3, 4}}, + {-2, {-3, 4}}, + {-1, {-3, 4}}, + {0, {-3, 4}}, + {1, {-3, 4}}, + {2, {-3, 4}}, + {3, {-3, 4}}, + {4, {-3, 4}}, + {5, {-3, 4}}, + {6, {-3, 4}}, + {7, {-3, 4}} + }; + + bool is_set_empty = true; + std::size_t ie = 0; + set( + [&expected, &is_set_empty, &ie](const auto& x_interval, const auto& yz) + { + is_set_empty = false; + EXPECT_EQ(expected[ie++], std::make_pair(yz[0], x_interval)); + }); + EXPECT_EQ(ie, expected.size()); + EXPECT_FALSE(is_set_empty); + } + { + const auto translated_ca = translate(ca, {3 + 1, 3 + 1}); + const auto joined_cas = union_(ca, translated_ca); + + const auto set = expand(joined_cas, 3); + + expected_t expected{ + {-3, {-3, 4}}, + {-2, {-3, 4}}, + {-1, {-3, 4}}, + {0, {-3, 4}}, + {1, {-3, 8}}, + {2, {-3, 8}}, + {3, {-3, 8}}, + {4, {1, 8} }, + {5, {1, 8} }, + {6, {1, 8} }, + {7, {1, 8} } + }; + + bool is_set_empty = true; + std::size_t ie = 0; + set( + [&expected, &is_set_empty, &ie](const auto& x_interval, const auto& yz) + { + is_set_empty = false; + EXPECT_EQ(expected[ie++], std::make_pair(yz[0], x_interval)); + }); + EXPECT_EQ(ie, expected.size()); + EXPECT_FALSE(is_set_empty); + + const auto lca_joined_cas = joined_cas.to_lca(); + const auto lca_set = set.to_lca(); + } + } + + TEST(subset, contract) + { + LevelCellArray<2> ca; + + ca.add_interval_back({0, 1}, {0}); + + { + const auto translated_ca = translate(ca, {3 + 1, 0}); + const auto joined_cas = union_(ca, translated_ca); + + const auto set = contract(joined_cas, 1); + + bool is_set_empty = true; + set( + [&is_set_empty](const auto& x_interval, const auto& yz) + { + fmt::print("x_interval = {} -- yz = {}", x_interval, yz[0]); + is_set_empty = false; + }); + EXPECT_TRUE(is_set_empty); + //~ EXPECT_TRUE(set.empty()); + } + } + TEST(subset, translate) { CellList<1> cl; @@ -1172,6 +1305,10 @@ namespace samurai ca = {cl, true}; + std::cout << self(ca[4]).on(3).to_lca() << std::endl; + + fmt::print("===============================================\n"); + // Test self-similarity at different scales bool found = false; apply(intersection(ca[3], self(ca[4]).on(3)),