tpl_r_tree.H

Go to the documentation of this file.

/*
                          Aleph_w
 
  Data structures & Algorithms
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/
 
#ifndef TPL_R_TREE_H
#define TPL_R_TREE_H
 
#include <algorithm>
#include <array>
#include <concepts>
#include <cstddef>
#include <limits>
#include <memory>
#include <type_traits>
#include <utility>
 
#include <ah-errors.H>
#include <point.H>
#include <tpl_array.H>
 
namespace Aleph {
 
enum class RTreeVariant
{
  Guttman,
  RStar
};
 
template <typename Payload, size_t MaxEntries = 16, size_t MinEntries = MaxEntries / 2, RTreeVariant Variant = RTreeVariant::Guttman>
class RTree
{
  static_assert(MaxEntries >= 2, "RTree requires MaxEntries >= 2");
  static_assert(MinEntries >= 1, "RTree requires MinEntries >= 1");
  static_assert(2 * MinEntries <= MaxEntries + 1, "RTree requires 2 * MinEntries <= MaxEntries + 1");
  static_assert(std::default_initializable<Payload>,
                "RTree requires default-initializable Payload");
  static_assert(std::movable<Payload>, "RTree requires movable Payload");
 
public:
  struct Entry
  {
    Rectangle bbox;  
    Payload value;   
  };
 
  struct DebugNode
  {
    Rectangle bbox;          
    bool is_leaf = false;    
    size_t depth = 0;        
    Array<size_t> children = Array<size_t>(0);  
    Array<Rectangle> entry_boxes = Array<Rectangle>(0);  
  };
 
  struct DebugSnapshot
  {
    Array<DebugNode> nodes;             
    size_t root = std::numeric_limits<size_t>::max();  
  };
 
private:
  struct Node;
 
  struct Child
  {
    Rectangle bbox;               
    std::unique_ptr<Node> child;  
  };
 
  struct Node
  {
    bool leaf = true;
    Array<Entry> data = Array<Entry>(0);      
    Array<Child> children = Array<Child>(0);  
  };
 
  std::unique_ptr<Node> root_;
  size_t size_ = 0;    
  size_t height_ = 0;  
 
  // R*-only transient state (empty/false while no insertion is in progress).
  bool rstar_reinsert_available_ = false;  
  Array<Entry> rstar_reinsert_buffer_;     
 
  void clear_rstar_reinsert_state() noexcept
  {
    rstar_reinsert_available_ = false;
    rstar_reinsert_buffer_.clear();
  }
 
  // ---- geometry helpers -------------------------------------------------
 
  [[nodiscard]] static Rectangle union_bbox(const Rectangle &a, const Rectangle &b)
  {
    const Geom_Number xmin = a.get_xmin() < b.get_xmin() ? a.get_xmin() : b.get_xmin();
    const Geom_Number ymin = a.get_ymin() < b.get_ymin() ? a.get_ymin() : b.get_ymin();
    const Geom_Number xmax = a.get_xmax() > b.get_xmax() ? a.get_xmax() : b.get_xmax();
    const Geom_Number ymax = a.get_ymax() > b.get_ymax() ? a.get_ymax() : b.get_ymax();
    return {xmin, ymin, xmax, ymax};
  }
 
  [[nodiscard]] static Geom_Number enlargement(const Rectangle &base, const Rectangle &added)
  {
    return union_bbox(base, added).area() - base.area();
  }
 
  [[nodiscard]] static Geom_Number overlap_area(const Rectangle &a, const Rectangle &b)
  {
    const Geom_Number xlo = a.get_xmin() > b.get_xmin() ? a.get_xmin() : b.get_xmin();
    const Geom_Number xhi = a.get_xmax() < b.get_xmax() ? a.get_xmax() : b.get_xmax();
    if (xhi <= xlo)
      return Geom_Number(0);
    const Geom_Number ylo = a.get_ymin() > b.get_ymin() ? a.get_ymin() : b.get_ymin();
    const Geom_Number yhi = a.get_ymax() < b.get_ymax() ? a.get_ymax() : b.get_ymax();
    if (yhi <= ylo)
      return Geom_Number(0);
    return (xhi - xlo) * (yhi - ylo);
  }
 
  [[nodiscard]] static size_t entry_count(const Node &node)
  {
    return node.leaf ? node.data.size() : node.children.size();
  }
 
  [[nodiscard]] static Rectangle compute_mbr(const Node &node)
  {
    if (node.leaf)
      {
        Rectangle acc = node.data(0).bbox;
        for (size_t i = 1; i < node.data.size(); ++i)
          acc = union_bbox(acc, node.data(i).bbox);
        return acc;
      }
    Rectangle acc = node.children(0).bbox;
    for (size_t i = 1; i < node.children.size(); ++i)
      acc = union_bbox(acc, node.children(i).bbox);
    return acc;
  }
 
  [[nodiscard]] static size_t choose_subtree(const Node &node, const Rectangle &bbox)
  {
    if constexpr (Variant == RTreeVariant::RStar)
      if (node.children(0).child->leaf)
        return rstar_choose_overlap(node, bbox);
 
    size_t best = 0;
    Geom_Number best_enlarge = enlargement(node.children(0).bbox, bbox);
    Geom_Number best_area = node.children(0).bbox.area();
    for (size_t i = 1; i < node.children.size(); ++i)
      {
        const Geom_Number e = enlargement(node.children(i).bbox, bbox);
        const Geom_Number a = node.children(i).bbox.area();
        if (e < best_enlarge or (e == best_enlarge and a < best_area))
          {
            best = i;
            best_enlarge = e;
            best_area = a;
          }
      }
    return best;
  }
 
  [[nodiscard]] static size_t rstar_choose_overlap(const Node &node, const Rectangle &bbox)
  {
    const size_t k = node.children.size();
    size_t best = 0;
    bool first = true;
    Geom_Number best_ov = 0;
    Geom_Number best_enl = 0;
    Geom_Number best_area = 0;
    for (size_t i = 0; i < k; ++i)
      {
        const Rectangle grown = union_bbox(node.children(i).bbox, bbox);
        Geom_Number ov = 0;
        for (size_t j = 0; j < k; ++j)
          if (j != i)
            ov += overlap_area(grown, node.children(j).bbox) -
              overlap_area(node.children(i).bbox, node.children(j).bbox);
        const Geom_Number enl = grown.area() - node.children(i).bbox.area();
        const Geom_Number area = node.children(i).bbox.area();
        if (first or ov < best_ov or
            (ov == best_ov and (enl < best_enl or (enl == best_enl and area < best_area))))
          {
            best = i;
            best_ov = ov;
            best_enl = enl;
            best_area = area;
            first = false;
          }
      }
    return best;
  }
 
  // ---- quadratic split (Guttman) ----------------------------------------
 
  template <typename EntryT>
  [[nodiscard]] static Array<EntryT> quadratic_split(Array<EntryT> &entries)
  {
    const size_t n = entries.size();
    Array<EntryT> group1(n);
    Array<EntryT> group2(n);
 
    Array<EntryT> items = std::move(entries);
    entries = Array<EntryT>();
 
    // PickSeeds: the pair wasting the most area if grouped together.
    size_t seed1 = 0;
    size_t seed2 = 1;
    Geom_Number worst =
      union_bbox(items(0).bbox, items(1).bbox).area() - items(0).bbox.area() - items(1).bbox.area();
    for (size_t i = 0; i < n; ++i)
      for (size_t j = i + 1; j < n; ++j)
        {
          const Geom_Number d = union_bbox(items(i).bbox, items(j).bbox).area() -
            items(i).bbox.area() - items(j).bbox.area();
          if (d > worst)
            {
              worst = d;
              seed1 = i;
              seed2 = j;
            }
        }
 
    std::array<bool, MaxEntries + 1> assigned{};
    Rectangle mbr1 = items(seed1).bbox;
    Rectangle mbr2 = items(seed2).bbox;
    assigned[seed1] = true;
    assigned[seed2] = true;
    group1.append(std::move(items(seed1)));
    group2.append(std::move(items(seed2)));
    size_t remaining = n - 2;
 
    auto dump_into = [&](Array<EntryT> &group)
      {
        for (size_t k = 0; k < n; ++k)
          if (not assigned[k])
            {
              assigned[k] = true;
              group.append(std::move(items(k)));
            }
        remaining = 0;
      };
 
    while (remaining > 0)
      {
        // Force the rest into a group that would otherwise miss MinEntries.
        if (group1.size() + remaining == MinEntries)
          {
            dump_into(group1);
            break;
          }
        if (group2.size() + remaining == MinEntries)
          {
            dump_into(group2);
            break;
          }
 
        // PickNext: entry with the strongest preference for one group.
        size_t pick = n;
        Geom_Number best_diff = 0;
        Geom_Number pick_d1 = 0;
        Geom_Number pick_d2 = 0;
        for (size_t k = 0; k < n; ++k)
          {
            if (assigned[k])
              continue;
            const Geom_Number d1 = enlargement(mbr1, items(k).bbox);
            const Geom_Number d2 = enlargement(mbr2, items(k).bbox);
            const Geom_Number diff = d1 > d2 ? d1 - d2 : d2 - d1;
            if (pick == n or diff > best_diff)
              {
                best_diff = diff;
                pick = k;
                pick_d1 = d1;
                pick_d2 = d2;
              }
          }
 
        // Assign to the group needing less enlargement (ties: smaller area,
        // then fewer entries).
        bool to_group1;
        if (pick_d1 < pick_d2)
          to_group1 = true;
        else if (pick_d2 < pick_d1)
          to_group1 = false;
        else if (mbr1.area() != mbr2.area())
          to_group1 = mbr1.area() < mbr2.area();
        else
          to_group1 = group1.size() <= group2.size();
 
        assigned[pick] = true;
        if (to_group1)
          {
            mbr1 = union_bbox(mbr1, items(pick).bbox);
            group1.append(std::move(items(pick)));
          }
        else
          {
            mbr2 = union_bbox(mbr2, items(pick).bbox);
            group2.append(std::move(items(pick)));
          }
        --remaining;
      }
 
    entries = std::move(group1);
    return group2;
  }
 
  template <typename EntryT>
  [[nodiscard]] static Array<EntryT> rstar_split(Array<EntryT> &entries)
  {
    const size_t n = entries.size();
    Array<EntryT> group1(n);
    Array<EntryT> group2(n);
 
    Array<EntryT> items = std::move(entries);
    entries = Array<EntryT>();
 
    // prefix/suffix are indexed up to and including n (prefix[n], suffix
    // accessed down to suffix[0]), and n can reach MaxEntries + 1, so both
    // need MaxEntries + 2 slots -- one more than the MaxEntries + 1 used for
    // per-entry arrays like `order`.
    auto make_mbr_cache = [&items, n](const std::array<size_t, MaxEntries + 1> &order)
      {
        std::array<Rectangle, MaxEntries + 2> prefix{};
        std::array<Rectangle, MaxEntries + 2> suffix{};
        prefix[1] = items(order[0]).bbox;
        for (size_t k = 2; k <= n; ++k)
          prefix[k] = union_bbox(prefix[k - 1], items(order[k - 1]).bbox);
        suffix[n - 1] = items(order[n - 1]).bbox;
        for (size_t k = n - 1; k > 0; --k)
          suffix[k - 1] = union_bbox(items(order[k - 1]).bbox, suffix[k]);
        return std::pair{prefix, suffix};
      };
 
    auto make_order = [&items, n](const int axis, const int edge)
      {
        std::array<size_t, MaxEntries + 1> order{};
        for (size_t i = 0; i < n; ++i)
          order[i] = i;
        std::sort(order.begin(), order.begin() + n,
                  [&items, axis, edge](const size_t a, const size_t b)
          {
            const Rectangle &ra = items(a).bbox;
            const Rectangle &rb = items(b).bbox;
            const Geom_Number la = axis == 0 ? ra.get_xmin() : ra.get_ymin();
            const Geom_Number ua = axis == 0 ? ra.get_xmax() : ra.get_ymax();
            const Geom_Number lb = axis == 0 ? rb.get_xmin() : rb.get_ymin();
            const Geom_Number ub = axis == 0 ? rb.get_xmax() : rb.get_ymax();
            const Geom_Number pa = edge == 0 ? la : ua;
            const Geom_Number pb = edge == 0 ? lb : ub;
            if (pa != pb)
              return pa < pb;
            return (edge == 0 ? ua : la) < (edge == 0 ? ub : lb);
          });
        return order;
      };
 
    // ChooseSplitAxis: minimize the summed margin over both edge sortings and
    // all valid distributions.
    int chosen_axis = 0;
    Geom_Number best_margin_sum = 0;
    for (int axis = 0; axis < 2; ++axis)
      {
        Geom_Number margin_sum = 0;
        for (int edge = 0; edge < 2; ++edge)
          {
            const std::array<size_t, MaxEntries + 1> order = make_order(axis, edge);
            const auto [prefix, suffix] = make_mbr_cache(order);
            for (size_t sz = MinEntries; sz <= n - MinEntries; ++sz)
              margin_sum += prefix[sz].perimeter() + suffix[sz].perimeter();
          }
        if (axis == 0 or margin_sum < best_margin_sum)
          {
            best_margin_sum = margin_sum;
            chosen_axis = axis;
          }
      }
 
    // ChooseSplitIndex: on the chosen axis, minimize overlap area, then area.
    std::array<size_t, MaxEntries + 1> best_order{};
    size_t best_sz = MinEntries;
    bool first = true;
    Geom_Number best_overlap = 0;
    Geom_Number best_area = 0;
    for (int edge = 0; edge < 2; ++edge)
      {
        const std::array<size_t, MaxEntries + 1> order = make_order(chosen_axis, edge);
        const auto [prefix, suffix] = make_mbr_cache(order);
        for (size_t sz = MinEntries; sz <= n - MinEntries; ++sz)
          {
            const Rectangle &g1 = prefix[sz];
            const Rectangle &g2 = suffix[sz];
            const Geom_Number ov = overlap_area(g1, g2);
            const Geom_Number ar = g1.area() + g2.area();
            if (first or ov < best_overlap or (ov == best_overlap and ar < best_area))
              {
                first = false;
                best_overlap = ov;
                best_area = ar;
                best_order = order;
                best_sz = sz;
              }
          }
      }
 
    for (size_t k = 0; k < best_sz; ++k)
      group1.append(std::move(items(best_order[k])));
    for (size_t k = best_sz; k < n; ++k)
      group2.append(std::move(items(best_order[k])));
    entries = std::move(group1);
    return group2;
  }
 
  // ---- insertion --------------------------------------------------------
 
  template <typename EntryT>
  [[nodiscard]] static Array<EntryT> split_entries(Array<EntryT> &entries)
  {
    if constexpr (Variant == RTreeVariant::RStar)
      return rstar_split(entries);
    else
      return quadratic_split(entries);
  }
 
  void reinsert_farthest(Node &node)
  {
    const size_t n = node.data.size();
    const Point centre = compute_mbr(node).center();
 
    std::array<size_t, MaxEntries + 1> order{};
    for (size_t i = 0; i < n; ++i)
      order[i] = i;
    std::sort(order.begin(), order.begin() + n, [&node, &centre](const size_t a, const size_t b)
      {
        return node.data(a).bbox.center().distance_squared_to(centre) >
          node.data(b).bbox.center().distance_squared_to(centre);
      });
 
    size_t p = (3 * MaxEntries) / 10;
    if (p < 1)
      p = 1;
    if (p > MaxEntries + 1 - MinEntries)
      p = MaxEntries + 1 - MinEntries;
 
    Array<Entry> kept(n - p);
    rstar_reinsert_buffer_.reserve(rstar_reinsert_buffer_.size() + p);
    for (size_t k = p; k < n; ++k)
      kept.append(std::move(node.data(order[k])));
    for (size_t k = 0; k < p; ++k)
      rstar_reinsert_buffer_.append(std::move(node.data(order[k])));
    node.data = std::move(kept);
  }
 
  [[nodiscard]] std::unique_ptr<Node> insert_descend(Node &node, Entry entry)
  {
    if (node.leaf)
      {
        node.data.reserve(MaxEntries + 1);
        node.data.append(std::move(entry));
        if (node.data.size() <= MaxEntries)
          return nullptr;
 
        if constexpr (Variant == RTreeVariant::RStar)
          if (rstar_reinsert_available_ and &node != root_.get())
            {
              rstar_reinsert_available_ = false;
              reinsert_farthest(node);
              return nullptr;
            }
 
        Array<Entry> group2 = split_entries(node.data);
        auto sibling = std::make_unique<Node>();
        sibling->leaf = true;
        sibling->data = std::move(group2);
        return sibling;
      }
 
    const size_t idx = choose_subtree(node, entry.bbox);
    std::unique_ptr<Node> split = insert_descend(*node.children(idx).child, std::move(entry));
    node.children(idx).bbox = compute_mbr(*node.children(idx).child);
    if (split == nullptr)
      return nullptr;
 
    Rectangle split_mbr = compute_mbr(*split);
    node.children.reserve(MaxEntries + 1);
    node.children.append(Child{std::move(split_mbr), std::move(split)});
    if (node.children.size() <= MaxEntries)
      return nullptr;
 
    Array<Child> group2 = split_entries(node.children);
    auto sibling = std::make_unique<Node>();
    sibling->leaf = false;
    sibling->children = std::move(group2);
    return sibling;
  }
 
  void insert_one(Entry entry)
  {
    if (root_ == nullptr)
      {
        auto root = std::make_unique<Node>();
        root->leaf = true;
        root->data.reserve(1);
        root->data.append(std::move(entry));
        root_ = std::move(root);
        height_ = 1;
        return;
      }
 
    std::unique_ptr<Node> split = insert_descend(*root_, std::move(entry));
    if (split == nullptr)
      return;
 
    auto new_root = std::make_unique<Node>();
    new_root->leaf = false;
    new_root->children.reserve(2);
    new_root->children.append(Child{compute_mbr(*root_), std::move(root_)});
    new_root->children.append(Child{compute_mbr(*split), std::move(split)});
    root_ = std::move(new_root);
    ++height_;
  }
 
  void add_entry(Entry entry)
  {
    if constexpr (Variant == RTreeVariant::RStar)
      {
        rstar_reinsert_available_ = true;
        insert_one(std::move(entry));
        while (not rstar_reinsert_buffer_.is_empty())
          insert_one(rstar_reinsert_buffer_.remove_last());
        clear_rstar_reinsert_state();
      }
    else
      insert_one(std::move(entry));
  }
 
  // ---- deletion ---------------------------------------------------------
 
  static void collect_data_entries(Node &node, Array<Entry> &out)
  {
    if (node.leaf)
      {
        for (size_t i = 0; i < node.data.size(); ++i)
          out.append(std::move(node.data(i)));
        return;
      }
    for (size_t i = 0; i < node.children.size(); ++i)
      collect_data_entries(*node.children(i).child, out);
  }
 
  [[nodiscard]] static size_t data_entry_count(const Node &node)
  {
    if (node.leaf)
      return node.data.size();
    size_t total = 0;
    for (size_t i = 0; i < node.children.size(); ++i)
      {
        const size_t subtree_count = data_entry_count(*node.children(i).child);
        ah_overflow_error_if(total > std::numeric_limits<size_t>::max() - subtree_count)
          << "RTree: data entry count would overflow";
        total += subtree_count;
      }
    return total;
  }
 
  void erase_descend(Node &node, const Rectangle &bbox, const Payload &value, Array<Entry> &orphans,
                     bool &removed)
    requires std::equality_comparable<Payload>
  {
    if (node.leaf)
      {
        for (size_t i = 0; i < node.data.size(); ++i)
          if (node.data(i).bbox == bbox and node.data(i).value == value)
            {
              Array<Entry> kept(node.data.size() - 1);
              for (size_t j = 0; j < node.data.size(); ++j)
                if (j != i)
                  kept.append(std::move(node.data(j)));
              node.data = std::move(kept);
              removed = true;
              return;
            }
        return;
      }
 
    for (size_t i = 0; i < node.children.size(); ++i)
      {
        Child &c = node.children(i);
        if (not box_contains_box(c.bbox, bbox))
          continue;
 
        erase_descend(*c.child, bbox, value, orphans, removed);
        if (not removed)
          continue;
 
        if (entry_count(*c.child) < MinEntries)
          {
            const size_t orphan_count = data_entry_count(*c.child);
            ah_overflow_error_if(orphans.size() >
                                 std::numeric_limits<size_t>::max() - orphan_count)
              << "RTree: orphan collection would overflow";
            orphans.reserve(orphans.size() + orphan_count);
            Array<Child> kept(node.children.size() - 1);
            collect_data_entries(*c.child, orphans);
            for (size_t j = 0; j < node.children.size(); ++j)
              if (j != i)
                kept.append(std::move(node.children(j)));
            node.children = std::move(kept);
          }
        else
          c.bbox = compute_mbr(*c.child);
        return;
      }
  }
 
  [[nodiscard]] static bool box_contains_box(const Rectangle &outer, const Rectangle &inner)
  {
    return outer.get_xmin() <= inner.get_xmin() and outer.get_xmax() >= inner.get_xmax() and
      outer.get_ymin() <= inner.get_ymin() and outer.get_ymax() >= inner.get_ymax();
  }
 
  // ---- queries ----------------------------------------------------------
 
  template <typename F>
  void for_each_intersecting_rec(const Node &node, const Rectangle &rect, F &f) const
  {
    if (node.leaf)
      {
        for (size_t i = 0; i < node.data.size(); ++i)
          if (node.data(i).bbox.intersects(rect))
            f(node.data(i).bbox, node.data(i).value);
        return;
      }
    for (size_t i = 0; i < node.children.size(); ++i)
      if (node.children(i).bbox.intersects(rect))
        for_each_intersecting_rec(*node.children(i).child, rect, f);
  }
 
  template <typename F>
  void for_each_containing_rec(const Node &node, const Point &p, F &f) const
  {
    if (node.leaf)
      {
        for (size_t i = 0; i < node.data.size(); ++i)
          if (node.data(i).bbox.contains(p))
            f(node.data(i).bbox, node.data(i).value);
        return;
      }
    for (size_t i = 0; i < node.children.size(); ++i)
      if (node.children(i).bbox.contains(p))
        for_each_containing_rec(*node.children(i).child, p, f);
  }
 
  // ---- copy / verify ----------------------------------------------------
 
  [[nodiscard]] static std::unique_ptr<Node> clone_node(const Node &node)
    requires (std::is_copy_constructible_v<Payload> and std::movable<Payload>)
  {
    auto copy = std::make_unique<Node>();
    copy->leaf = node.leaf;
    if (node.leaf)
      {
        copy->data.reserve(node.data.size());
        for (size_t i = 0; i < node.data.size(); ++i)
          copy->data.append(Entry{node.data(i).bbox, node.data(i).value});
      }
    else
      {
        copy->children.reserve(node.children.size());
        for (size_t i = 0; i < node.children.size(); ++i)
          copy->children.append(Child{node.children(i).bbox, clone_node(*node.children(i).child)});
      }
    return copy;
  }
 
  [[nodiscard]] bool verify_rec(const Node &node, const size_t depth, const bool is_root,
                                size_t &leaf_depth, size_t &count, Rectangle &out_mbr) const
  {
    const size_t n = entry_count(node);
    if (n > MaxEntries)
      return false;
    if (not is_root and n < MinEntries)
      return false;
    if (n == 0)
      return false;  // an empty node must never be retained
 
    if (node.leaf)
      {
        if (leaf_depth == std::numeric_limits<size_t>::max())
          leaf_depth = depth;
        else if (leaf_depth != depth)
          return false;
        out_mbr = compute_mbr(node);
        if (count > std::numeric_limits<size_t>::max() - n)
          return false;
        count += n;
        return true;
      }
 
    if (is_root and n < 2)
      return false;  // a non-leaf root must have at least two children
 
    Rectangle acc;
    for (size_t i = 0; i < node.children.size(); ++i)
      {
        Rectangle child_mbr;
        if (not verify_rec(*node.children(i).child, depth + 1, false, leaf_depth, count, child_mbr))
          return false;
        if (node.children(i).bbox != child_mbr)
          return false;  // stored MBR must be tight
        acc = i == 0 ? child_mbr : union_bbox(acc, child_mbr);
      }
    out_mbr = acc;
    return true;
  }
 
public:
  RTree() noexcept = default;
 
  RTree(RTree &&other) noexcept
    : root_(std::move(other.root_)), size_(other.size_), height_(other.height_)
  {
    other.size_ = 0;
    other.height_ = 0;
    other.clear_rstar_reinsert_state();
  }
 
  RTree &operator = (RTree &&other) noexcept
  {
    if (this == &other)
      return *this;
    root_ = std::move(other.root_);
    size_ = other.size_;
    height_ = other.height_;
    clear_rstar_reinsert_state();
    other.size_ = 0;
    other.height_ = 0;
    other.clear_rstar_reinsert_state();
    return *this;
  }
 
  RTree(const RTree &other)
    requires (std::is_copy_constructible_v<Payload> and std::movable<Payload>)
    : size_(other.size_), height_(other.height_)
  {
    if (other.root_ != nullptr)
      root_ = clone_node(*other.root_);
  }
 
  RTree &operator = (const RTree &other)
    requires (std::is_copy_constructible_v<Payload> and std::movable<Payload>)
  {
    if (this == &other)
      return *this;
    std::unique_ptr<Node> new_root;
    if (other.root_ != nullptr)
      new_root = clone_node(*other.root_);
    root_ = std::move(new_root);
    size_ = other.size_;
    height_ = other.height_;
    return *this;
  }
 
  [[nodiscard]] bool is_empty() const noexcept
  {
    return size_ == 0;
  }
 
  [[nodiscard]] size_t size() const noexcept
  {
    return size_;
  }
 
  [[nodiscard]] size_t height() const noexcept
  {
    return height_;
  }
 
  void clear() noexcept
  {
    root_.reset();
    size_ = 0;
    height_ = 0;
    clear_rstar_reinsert_state();
  }
 
  void insert(const Rectangle &bbox, const Payload &value)
  {
    ah_overflow_error_if(size_ == std::numeric_limits<size_t>::max())
      << "RTree: size would overflow";
    add_entry(Entry{bbox, value});
    ++size_;
  }
 
  void insert(const Rectangle &bbox, Payload &&value)
  {
    ah_overflow_error_if(size_ == std::numeric_limits<size_t>::max())
      << "RTree: size would overflow";
    add_entry(Entry{bbox, std::move(value)});
    ++size_;
  }
 
  [[nodiscard]] bool erase(const Rectangle &bbox, const Payload &value)
    requires std::equality_comparable<Payload>
  {
    if (root_ == nullptr)
      return false;
 
    Array<Entry> orphans;
    bool removed = false;
    erase_descend(*root_, bbox, value, orphans, removed);
    if (not removed)
      return false;
 
    --size_;
 
    // Reinsert the entries orphaned by condensing (they remain in the tree, so
    // add_entry must not change size_).
    for (size_t i = 0; i < orphans.size(); ++i)
      add_entry(std::move(orphans(i)));
 
    // Shrink a non-leaf root that ended up with a single child.
    while (root_ != nullptr and not root_->leaf and root_->children.size() == 1)
      {
        root_ = std::move(root_->children(0).child);
        --height_;
      }
    // Drop an emptied leaf root.
    if (root_ != nullptr and root_->leaf and root_->data.is_empty())
      {
        root_.reset();
        height_ = 0;
      }
    return true;
  }
 
  template <typename F>
  void for_each_intersecting(const Rectangle &rect, F &&f) const
  {
    if (root_ != nullptr)
      for_each_intersecting_rec(*root_, rect, f);
  }
 
  [[nodiscard]] Array<Payload> search_intersects(const Rectangle &rect) const
    requires (std::is_copy_constructible_v<Payload> and std::movable<Payload>)
  {
    Array<Payload> out;
    for_each_intersecting(rect, [&out](const Rectangle &, const Payload &v)
      {
        out.append(Payload(v));
      });
    return out;
  }
 
  [[nodiscard]] Array<Payload> search_contains(const Point &p) const
    requires (std::is_copy_constructible_v<Payload> and std::movable<Payload>)
  {
    Array<Payload> out;
    auto collect = [&out](const Rectangle &, const Payload &v)
      {
        out.append(Payload(v));
      };
    if (root_ != nullptr)
      for_each_containing_rec(*root_, p, collect);
    return out;
  }
 
  [[nodiscard]] bool verify() const
  {
    if (root_ == nullptr)
      return size_ == 0 and height_ == 0;
 
    size_t leaf_depth = std::numeric_limits<size_t>::max();
    size_t count = 0;
    Rectangle mbr;
    if (not verify_rec(*root_, 1, true, leaf_depth, count, mbr))
      return false;
    return count == size_ and leaf_depth == height_;
  }
 
  [[nodiscard]] DebugSnapshot debug_snapshot() const
  {
    DebugSnapshot snap;
    if (root_ == nullptr)
      return snap;
 
    // Each `snap.nodes(out_idx)` access below is by index, re-fetched after
    // every recursive call: `Array::append` may reallocate the backing
    // storage, so a reference taken before recursing could dangle.
    auto dfs = [&](const auto &self, const Node &node, const size_t depth) -> size_t
      {
        const size_t out_idx = snap.nodes.size();
        snap.nodes.append(DebugNode{});
        snap.nodes(out_idx).is_leaf = node.leaf;
        snap.nodes(out_idx).depth = depth;
 
        if (node.leaf)
          {
            // Safe to populate in place: this loop does not recurse, so
            // `snap.nodes` cannot reallocate between the reference below
            // and its last use.
            Array<Rectangle> &boxes = snap.nodes(out_idx).entry_boxes;
            boxes.reserve(node.data.size());
            for (size_t i = 0; i < node.data.size(); ++i)
              boxes.append(node.data(i).bbox);
          }
        else
          {
            snap.nodes(out_idx).children.reserve(node.children.size());
            for (size_t i = 0; i < node.children.size(); ++i)
              {
                const size_t child_idx = self(self, *node.children(i).child, depth + 1);
                snap.nodes(out_idx).children.append(child_idx);
              }
          }
 
        snap.nodes(out_idx).bbox = compute_mbr(node);
        return out_idx;
      };
 
    snap.root = dfs(dfs, *root_, 0);
    return snap;
  }
};
 
}  // namespace Aleph
 
#endif  // TPL_R_TREE_H