Commit f804ed18 authored by Edwin Carlinet's avatar Edwin Carlinet
Browse files

Faster immersion using hvectors for hierarchical queues.

parent 490c79ad
......@@ -65,6 +65,7 @@ target_sources(Pylene PRIVATE
src/morpho/immersion.cpp
src/morpho/component_tree.cpp
src/morpho/hvector.cpp
src/morpho/hvector_unbounded.cpp
)
# Compiler configurations
......
......@@ -137,6 +137,8 @@ namespace mln::morpho::experimental
std::pair<component_tree<image_value_t<I>>, image_ch_value_t<I, component_tree<>::node_id_type>> //
maxtree(I input, N nbh)
{
mln_entering("mln::morpho::experimental::maxtree");
using V = image_value_t<I>;
constexpr std::size_t kStackReserve = 256;
......
#pragma once
#include <mln/core/concept/new/images.hpp>
#include <mln/core/assert.hpp>
#include <mln/core/algorithm/fill.hpp>
namespace mln::morpho::experimental::details
{
/// \brief Set of points ordered hierarchically into a simple list of points (only access to the front)
///
/// \tparam P Type of point
/// \tparam LinkImage Type of the link image
template <int N, class P, class LinkImage>
class hsimple_linked_lists
{
static_assert(N <= (1 << 24), "Too many numbers of level");
static_assert(std::is_convertible_v<P, image_point_t<LinkImage>>);
template <class J>
hsimple_linked_lists(J&& f);
template <class J>
hsimple_linked_lists(J&& f, P sentinel);
void push_front(int level, P p) noexcept;
P pop_front(int level) noexcept;
P front(int level) const noexcept;
bool empty(int level) const noexcept;
/// Get the level 𝑙 of the lowest non-empty queue with level ≤ 𝑙
/// such that every queue at level x with level ≤ x < 𝑙 are empty.
int lower_bound(int level) const noexcept;
/// Get the level 𝑙 of the highest non-empty queue with 𝑙 < level
/// such that every queue at level x with 𝑙 < x ≤ level are empty.
int upper_bound(int level) const noexcept;
private:
const P m_sentinel;
std::array<P, N> m_heads;
LinkImage m_next;
};
/******************************************/
/**** Implementation ****/
/******************************************/
template <int N, class P, class LinkImage>
template <class J>
inline hsimple_linked_lists<N, P, LinkImage>::hsimple_linked_lists(J&& f, P sentinel)
: m_sentinel{sentinel}
, m_next(std::forward<J>(f), image_build_params{})
{
m_heads.fill(m_sentinel);
}
template <int N, class P, class LinkImage>
template <class J>
inline hsimple_linked_lists<N, P, LinkImage>::hsimple_linked_lists(J&& f)
: hsimple_linked_lists{std::forward<J>(f), f.domain().br()}
{
}
template <int N, class P, class LinkImage>
inline bool hsimple_linked_lists<N, P, LinkImage>::empty(int level) const noexcept
{
mln_precondition(level < N);
return m_heads[level] == m_sentinel;
}
template <int N, class P, class LinkImage>
inline void hsimple_linked_lists<N, P, LinkImage>::push_front(int level, P p) noexcept
{
mln_precondition(level < N);
m_next(p) = m_heads[level];
m_heads[level] = p;
}
template <int N, class P, class LinkImage>
inline P hsimple_linked_lists<N, P, LinkImage>::pop_front(int level) noexcept
{
mln_precondition(level < N);
mln_precondition(m_heads[level] != m_sentinel && "Empty list");
P head = m_heads[level];
m_heads[level] = m_next(head);
return head;
}
template <int N, class P, class LinkImage>
inline P hsimple_linked_lists<N, P, LinkImage>::front(int level) const noexcept
{
mln_precondition(level < N);
mln_precondition(m_heads[level].size > 0 && "Empty list");
return m_heads[level].head;
}
template <int N, class P, class LinkImage>
inline int hsimple_linked_lists<N, P, LinkImage>::lower_bound(int level) const noexcept
{
mln_precondition(level < N);
while (level < N && m_heads[level] == m_sentinel)
++level;
return level;
}
template <int N, class P, class LinkImage>
inline int hsimple_linked_lists<N, P, LinkImage>::upper_bound(int level) const noexcept
{
mln_precondition(level < N);
while (level >= 0 && m_heads[level] == m_sentinel)
--level;
return level;
}
} // namespace mln::morpho::experimental::detail
#pragma once
#include <mln/core/assert.hpp>
#include <mln/core/concept/new/images.hpp>
#include <mln/data/experimental/histogram.hpp>
#include <algorithm>
#include <numeric>
namespace mln::morpho::experimental::details
{
/// \brief Set of points ordered hierarchically into a list of points (vector implemented)
///
/// \tparam P Type of point
template <class P>
class hvectors_unbounded;
/// \brief Base implementation
template <>
class hvectors_unbounded<void>
{
public:
bool empty(int level) const noexcept;
/// Get the level 𝑙 of the lowest non-empty queue with level ≤ 𝑙
/// such that every queue at level x with level ≤ x < 𝑙 are empty.
int lower_bound(int level) const noexcept;
/// Get the level 𝑙 of the highest non-empty queue with 𝑙 < level
/// such that every queue at level x with 𝑙 < x ≤ level are empty.
int upper_bound(int level) const noexcept;
protected:
/// \param nlevel Number of levels
/// \param size Size in bytes of an element
hvectors_unbounded(int nlevel, std::size_t size);
/// Realloc the buffer at a given level to double the capacity
/// \param level The level to free
/// \param size Size in bytes of an element
void resize(int level, std::size_t size);
/// Free the memory
void release();
virtual ~hvectors_unbounded() = default;
virtual void uninitialized_copy_n(void* in, std::size_t count, void* out) = 0;
virtual void destroy_n(void* buffer, std::size_t count) = 0;
struct node_t
{
void* begin;
int size;
int capacity;
};
node_t* m_lists; /* List of unitialized storage */
int m_nlevels;
};
template <class P>
class hvectors_unbounded : public hvectors_unbounded<void>
{
hvectors_unbounded(int nlevel);
~hvectors_unbounded() final;
hvectors_unbounded(const hvectors_unbounded&) = delete;
hvectors_unbounded(hvectors_unbounded&&) = delete;
hvectors_unbounded& operator= (const hvectors_unbounded&) = delete;
hvectors_unbounded& operator= (hvectors_unbounded&&) = delete;
void push_front(int level, P p) noexcept;
P pop_front(int level) noexcept;
P back(int level) const noexcept;
P front(int level) const noexcept;
private:
void uninitialized_copy_n(void* in, std::size_t count, void* out) final;
void destroy_n(void* buffer, std::size_t count) final;
using base = hvectors_unbounded<void>;
};
/******************************************/
/**** Implementation ****/
/******************************************/
inline int hvectors_unbounded<void>::lower_bound(int level) const noexcept
{
mln_precondition(level < m_nlevels);
while (level < m_nlevels && m_lists[level].size == 0)
++level;
return level;
}
inline int hvectors_unbounded<void>::upper_bound(int level) const noexcept
{
mln_precondition(level < m_nlevels);
while (level >= 0 && m_lists[level].size == 0)
--level;
return level;
}
inline bool hvectors_unbounded<void>::empty(int level) const noexcept
{
return m_lists[level].size == 0;
}
template <class P>
inline hvectors_unbounded<P>::hvectors_unbounded(int nlevel)
: base{nlevel, sizeof(P)}
{
}
template <class P>
inline hvectors_unbounded<P>::~hvectors_unbounded()
{
this->release();
}
template <class P>
void hvectors_unbounded<P>::uninitialized_copy_n(void* in, std::size_t count, void* out)
{
std::uninitialized_copy_n((P*)in, count, (P*)out);
}
template <class P>
void hvectors_unbounded<P>::destroy_n(void* in, std::size_t count)
{
std::destroy_n((P*)in, count);
}
template <class P>
inline void hvectors_unbounded<P>::push_front(int level, P p) noexcept
{
if (m_lists[level].size == m_lists[level].capacity)
this->resize(level, sizeof(P));
P* buffer = reinterpret_cast<P*>(m_lists[level].begin);
P* e = buffer + m_lists[level].size++;
new ((void*)e) P(p);
}
template <class P>
inline P hvectors_unbounded<P>::pop_front(int level) noexcept
{
mln_precondition(!this->empty(level) && "Empty list");
P* buffer = reinterpret_cast<P*>(m_lists[level].begin);
--m_lists[level].size;
P tmp = std::move(buffer[m_lists[level].size]);
std::destroy_at(buffer + m_lists[level].size);
return tmp;
}
template <class P>
inline P hvectors_unbounded<P>::front(int level) const noexcept
{
mln_precondition(!this->empty(level) && "Empty list");
P* buffer = reinterpret_cast<P*>(m_lists[level].begin);
return buffer[m_lists[level].size - 1];
}
} // namespace mln::morpho::experimental::detail
......@@ -2,6 +2,7 @@
#include <mln/core/concept/new/images.hpp>
#include <mln/core/algorithm/fill.hpp>
#include <mln/core/algorithm/transform.hpp>
#include <mln/core/extension/border_management.hpp>
#include <mln/core/extension/fill.hpp>
#include <mln/core/trace.hpp>
......@@ -22,6 +23,25 @@ namespace mln::morpho::experimental::details
/**** Implementation ****/
/******************************************/
template <class V>
struct irange
{
V inf;
V sup;
};
template <class I>
[[gnu::noinline]]
auto to_infsup(I inf, I sup)
{
mln_entering("mln::morpho::details::to_infsup");
using V = image_value_t<I>;
image_ch_value_t<I, irange<V>> out = imchvalue<irange<V>>(inf);
mln::transform(inf, sup, out, [](V a, V b) -> irange<V> { return {a, b}; });
return out;
}
template <class I>
void propagation(I inf, I sup, image_ch_value_t<I, int> ord, image_point_t<I> pstart, int& max_depth)
......@@ -58,12 +78,13 @@ namespace mln::morpho::experimental::details
// if (compute_indexes)
// sorted_indexes->reserve(ord.domain().size());
auto F = to_infsup(inf, sup);
pset<I> queue(inf);
auto p = pstart;
V previous_level = inf(p);
queue.insert(previous_level, p);
ord(p) = 0;
// if (compute_indexes)
// sorted_indexes->push_back(p);
......@@ -71,40 +92,50 @@ namespace mln::morpho::experimental::details
int depth = 0;
while (!queue.empty())
{
V cur_level;
P p;
auto tmp = queue.try_pop(previous_level);
if (tmp.has_value())
{
cur_level = previous_level;
p = tmp.value();
}
else
{
std::tie(cur_level, p) = queue.pop(previous_level);
auto [cur_level, p] = queue.pop(previous_level);
if (cur_level != previous_level)
++depth;
}
ord(p) = depth;
// if (compute_indexes)
// sorted_indexes->push_back(p);
for (auto q : nbh(p))
bool has_nbh_at_level;
while (true)
{
if (ord.at(q) == UNPROCESSED)
P nextp;
has_nbh_at_level = false;
ord(p) = depth;
for (auto q : nbh(p))
{
auto m = inf(q);
auto M = sup(q);
if (ord.at(q) != UNPROCESSED)
continue;
auto [m, M] = F(q);
if (M < cur_level)
queue.insert(M, q);
else if (cur_level < m)
queue.insert(m, q);
else
else if (has_nbh_at_level)
queue.insert(cur_level, q);
else
{
has_nbh_at_level = true;
nextp = q;
continue;
}
ord(q) = PROCESSED;
}
// If propagation stops at this point
if (!has_nbh_at_level)
break;
p = nextp;
}
previous_level = cur_level;
}
max_depth = depth;
......
#pragma once
#include <mln/morpho/experimental/private/hlinked_lists.hpp>
#include <mln/morpho/experimental/private/hsimple_linked_lists.hpp>
#include <mln/morpho/experimental/private/hvector_unbounded.hpp>
#include <limits>
#include <optional>
......@@ -44,8 +45,9 @@ namespace mln::morpho::experimental::details
static_assert(std::numeric_limits<Key>::digits <= 16);
static constexpr int nLevel = 1 << std::numeric_limits<Key>::digits;
using link_image_t = image_ch_value_t<I, Point>;
using delegate_t = hlinked_lists<nLevel, Point, link_image_t>;
//using link_image_t = image_ch_value_t<I, Point>;
//using delegate_t = hsimple_linked_lists<nLevel, Point, link_image_t>;
using delegate_t = hvectors_unbounded<Point>;
delegate_t m_delegate;
std::size_t m_size;
......@@ -58,8 +60,8 @@ namespace mln::morpho::experimental::details
template <class I>
template <class J>
pset<I>::pset(J&& ima)
: m_delegate{std::forward<J>(ima)},
pset<I>::pset(J&&)
: m_delegate{nLevel},
m_size{0}
{
}
......
#include <mln/morpho/experimental/private/hvector_unbounded.hpp>
namespace mln::morpho::experimental::details
{
namespace
{
static constexpr int kInitReservationPerLevel = 256;
}
hvectors_unbounded<void>::hvectors_unbounded(int nlevel, std::size_t size)
{
constexpr std::size_t K = alignof(std::max_align_t);
std::size_t hdr_size = sizeof(node_t) * nlevel;
std::size_t data_size = size * nlevel * kInitReservationPerLevel;
std::size_t hdr_full_size = (hdr_size + (K - 1)) / K * K; // header with padding for alignment
std::size_t total_size = hdr_full_size + data_size;
assert(hdr_full_size % K == 0 && hdr_size <= hdr_full_size && hdr_full_size < hdr_size + K);
m_nlevels = nlevel;
std::byte* buffer = (std::byte*)std::malloc(total_size);
m_lists = reinterpret_cast<node_t*>(buffer);
for (int i = 0; i < nlevel; ++i)
{
m_lists[i].size = 0;
m_lists[i].capacity = kInitReservationPerLevel;
m_lists[i].begin = (void*) (buffer + hdr_full_size + size * kInitReservationPerLevel * i);
}
}
void hvectors_unbounded<void>::release()
{
for (int i = 0; i < m_nlevels; ++i)
{
this->destroy_n(m_lists[i].begin, m_lists[i].size);
if (m_lists[i].capacity > kInitReservationPerLevel)
std::free(m_lists[i].begin);
}
std::free(m_lists);
}
void hvectors_unbounded<void>::resize(int i, std::size_t size)
{
constexpr float kResizingFactor = 2;
int old_capacity = m_lists[i].capacity;
void* old_ptr = m_lists[i].begin;
m_lists[i].capacity = (int)(old_capacity * kResizingFactor);
m_lists[i].begin = std::malloc(size * m_lists[i].capacity);
this->uninitialized_copy_n(old_ptr, m_lists[i].size, m_lists[i].begin);
this->destroy_n(old_ptr, m_lists[i].size);
if (old_capacity > kInitReservationPerLevel)
std::free(old_ptr);
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment