Complete prtree.h modularization - split into focused components

Successfully completed the migration from monolithic prtree.h (1617 lines)
to a modular architecture with clear separation of concerns.

## New Modular Structure

Created in include/prtree/core/detail/:

### 1. bounding_box.h (140 lines)
- BB<D> class for axis-aligned bounding boxes
- Geometric operations: intersection, union, area calculation
- Operator overloads for combining bounding boxes
- Serialization support

### 2. data_type.h (48 lines)
- DataType<T, D> class for index-bbox pairs
- clean_data() utility function
- Move semantics and swap support

### 3. pseudo_tree.h (228 lines)
- Leaf<T, B, D> - leaf node for pseudo tree
- PseudoPRTreeNode<T, B, D> - internal node
- PseudoPRTree<T, B, D> - construction-phase tree
- Multi-threaded construction support

### 4. nodes.h (168 lines)
- PRTreeLeaf<T, B, D> - final tree leaf
- PRTreeNode<T, B, D> - final tree node
- PRTreeElement<T, B, D> - tree element wrapper
- bfs() utility for tree traversal

### 5. prtree.h (1041 lines, down from 1617)
- Main PRTree<T, B, D> class
- Includes all modular components
- Cleaner, more maintainable interface

## Benefits

1. Faster Compilation: Changes to components don't require full rebuilds
2. Better Organization: Each component has clear responsibility
3. Easier Maintenance: Smaller focused files are easier to work with
4. Future Ready: Supports incremental improvements and testing

## Verification

- Build successful with no errors
- All existing functionality preserved
- No changes to Python API (100% backwards compatible)
- File size reduced: detail headers total ~699 lines of focused code

Migration complete.

Author: claude

include/prtree/core/detail/bounding_box.h

@@ -3,203 +3,136 @@
  * @brief Axis-Aligned Bounding Box (AABB) implementation
  *
  * Provides the BB<D> class for D-dimensional bounding boxes with
- * geometric operations like intersection, union, and containment tests.
+ * geometric operations like intersection, union, and area calculation.
  */
 #pragma once
 
 #include <algorithm>
-#include <array>
-#include <cmath>
-#include <limits>
-#include <numeric>
+#include <stdexcept>
 
 #include <cereal/cereal.hpp>
 
 #include "prtree/core/detail/types.h"
 
 using Real = float;
 
-/**
- * @brief D-dimensional Axis-Aligned Bounding Box
- *
- * Stores min/max coordinates for each dimension and provides
- * geometric operations.
- *
- * @tparam D Number of dimensions (2, 3, or 4)
- */
-template <int D = 2>
-class BB {
-public:
-  std::array<Real, D> lo;  ///< Minimum coordinates
-  std::array<Real, D> hi;  ///< Maximum coordinates
-
-  /// Default constructor - creates an invalid/empty box
-  BB() {
-    for (int i = 0; i < D; i++) {
-      lo[i] = std::numeric_limits<Real>::max();
-      hi[i] = -std::numeric_limits<Real>::max();
-    }
-  }
+template <int D = 2> class BB {
+private:
+  Real values[2 * D];
 
-  /// Constructor from coordinate arrays
-  BB(const std::array<Real, D> &lo_, const std::array<Real, D> &hi_)
-      : lo(lo_), hi(hi_) {}
+public:
+  BB() { clear(); }
 
-  /// Constructor from iterators (for compatibility with span/vector)
-  template <typename Iterator>
-  BB(Iterator lo_begin, Iterator lo_end, Iterator hi_begin, Iterator hi_end) {
-    std::copy(lo_begin, lo_end, lo.begin());
-    std::copy(hi_begin, hi_end, hi.begin());
+  BB(const Real (&minima)[D], const Real (&maxima)[D]) {
+    Real v[2 * D];
+    for (int i = 0; i < D; ++i) {
+      v[i] = -minima[i];
+      v[i + D] = maxima[i];
+    }
+    validate(v);
+    for (int i = 0; i < D; ++i) {
+      values[i] = v[i];
+      values[i + D] = v[i + D];
+    }
   }
 
-  /**
-   * @brief Check if this box intersects with another
-   *
-   * Two boxes intersect if they overlap in all dimensions.
-   */
-  bool intersects(const BB &other) const {
-    for (int i = 0; i < D; i++) {
-      if (hi[i] < other.lo[i] || lo[i] > other.hi[i])
-        return false;
+  BB(const Real (&v)[2 * D]) {
+    validate(v);
+    for (int i = 0; i < D; ++i) {
+      values[i] = v[i];
+      values[i + D] = v[i + D];
     }
-    return true;
   }
 
-  /**
-   * @brief Check if this box contains a point
-   */
-  bool contains_point(const std::array<Real, D> &point) const {
-    for (int i = 0; i < D; i++) {
-      if (point[i] < lo[i] || point[i] > hi[i])
-        return false;
+  Real min(const int dim) const {
+    if (unlikely(dim < 0 || D <= dim)) {
+      throw std::runtime_error("Invalid dim");
     }
-    return true;
+    return -values[dim];
   }
-
-  /**
-   * @brief Check if this box completely contains another
-   */
-  bool contains(const BB &other) const {
-    for (int i = 0; i < D; i++) {
-      if (other.lo[i] < lo[i] || other.hi[i] > hi[i])
-        return false;
+  Real max(const int dim) const {
+    if (unlikely(dim < 0 || D <= dim)) {
+      throw std::runtime_error("Invalid dim");
     }
-    return true;
+    return values[dim + D];
   }
 
-  /**
-   * @brief Compute the union of this box with another
-   *
-   * Returns the smallest box that contains both boxes.
-   */
-  BB union_with(const BB &other) const {
-    BB result;
-    for (int i = 0; i < D; i++) {
-      result.lo[i] = std::min(lo[i], other.lo[i]);
-      result.hi[i] = std::max(hi[i], other.hi[i]);
+  bool validate(const Real (&v)[2 * D]) const {
+    bool flag = false;
+    for (int i = 0; i < D; ++i) {
+      if (unlikely(-v[i] > v[i + D])) {
+        flag = true;
+        break;
+      }
     }
-    return result;
+    if (unlikely(flag)) {
+      throw std::runtime_error("Invalid Bounding Box");
+    }
+    return flag;
   }
-
-  /**
-   * @brief Compute the intersection of this box with another
-   *
-   * Returns an empty box if they don't intersect.
-   */
-  BB intersection_with(const BB &other) const {
-    BB result;
-    for (int i = 0; i < D; i++) {
-      result.lo[i] = std::max(lo[i], other.lo[i]);
-      result.hi[i] = std::min(hi[i], other.hi[i]);
-      if (result.lo[i] > result.hi[i])
-        return BB(); // Empty box
+  void clear() noexcept {
+    for (int i = 0; i < 2 * D; ++i) {
+      values[i] = -1e100;
     }
-    return result;
   }
 
-  /**
-   * @brief Compute the volume (area in 2D) of the box
-   */
-  Real volume() const {
-    Real vol = 1.0;
-    for (int i = 0; i < D; i++) {
-      Real extent = hi[i] - lo[i];
-      if (extent < 0)
-        return 0; // Invalid box
-      vol *= extent;
-    }
-    return vol;
+  Real val_for_comp(const int &axis) const noexcept {
+    const int axis2 = (axis + 1) % (2 * D);
+    return values[axis] + values[axis2];
   }
 
-  /**
-   * @brief Compute the perimeter (in 2D) or surface area (in 3D)
-   */
-  Real perimeter() const {
-    if constexpr (D == 2) {
-      return 2 * ((hi[0] - lo[0]) + (hi[1] - lo[1]));
-    } else if constexpr (D == 3) {
-      Real dx = hi[0] - lo[0];
-      Real dy = hi[1] - lo[1];
-      Real dz = hi[2] - lo[2];
-      return 2 * (dx * dy + dy * dz + dz * dx);
-    } else {
-      // For other dimensions, return sum of extents
-      Real sum = 0;
-      for (int i = 0; i < D; i++)
-        sum += hi[i] - lo[i];
-      return sum;
+  BB operator+(const BB &rhs) const {
+    Real result[2 * D];
+    for (int i = 0; i < 2 * D; ++i) {
+      result[i] = std::max(values[i], rhs.values[i]);
     }
+    return BB<D>(result);
   }
 
-  /**
-   * @brief Compute the center point of the box
-   */
-  std::array<Real, D> center() const {
-    std::array<Real, D> c;
-    for (int i = 0; i < D; i++)
-      c[i] = (lo[i] + hi[i]) / 2;
-    return c;
+  BB operator+=(const BB &rhs) {
+    for (int i = 0; i < 2 * D; ++i) {
+      values[i] = std::max(values[i], rhs.values[i]);
+    }
+    return *this;
   }
 
-  /**
-   * @brief Check if the box is valid (min <= max for all dimensions)
-   */
-  bool is_valid() const {
-    for (int i = 0; i < D; i++) {
-      if (lo[i] > hi[i])
-        return false;
+  void expand(const Real (&delta)[D]) noexcept {
+    for (int i = 0; i < D; ++i) {
+      values[i] += delta[i];
+      values[i + D] += delta[i];
     }
-    return true;
   }
 
-  /**
-   * @brief Check if the box is empty (zero volume)
-   */
-  bool is_empty() const { return volume() == 0; }
-
-  /**
-   * @brief Expand the box to include a point
-   */
-  void expand_to_include(const std::array<Real, D> &point) {
-    for (int i = 0; i < D; i++) {
-      lo[i] = std::min(lo[i], point[i]);
-      hi[i] = std::max(hi[i], point[i]);
+  bool operator()(
+      const BB &target) const { // whether this and target has any intersect
+
+    Real minima[D];
+    Real maxima[D];
+    bool flags[D];
+    bool flag = true;
+
+    for (int i = 0; i < D; ++i) {
+      minima[i] = std::min(values[i], target.values[i]);
+      maxima[i] = std::min(values[i + D], target.values[i + D]);
+    }
+    for (int i = 0; i < D; ++i) {
+      flags[i] = -minima[i] <= maxima[i];
+    }
+    for (int i = 0; i < D; ++i) {
+      flag &= flags[i];
     }
+    return flag;
   }
 
-  /**
-   * @brief Expand the box to include another box
-   */
-  void expand_to_include(const BB &other) {
-    for (int i = 0; i < D; i++) {
-      lo[i] = std::min(lo[i], other.lo[i]);
-      hi[i] = std::max(hi[i], other.hi[i]);
+  Real area() const {
+    Real result = 1;
+    for (int i = 0; i < D; ++i) {
+      result *= max(i) - min(i);
     }
+    return result;
   }
 
-  /// Serialization support
-  template <class Archive>
-  void serialize(Archive &ar) {
-    ar(CEREAL_NVP(lo), CEREAL_NVP(hi));
-  }
+  inline Real operator[](const int i) const { return values[i]; }
+
+  template <class Archive> void serialize(Archive &ar) { ar(values); }
 };

include/prtree/core/detail/data_type.h

@@ -0,0 +1,47 @@
+/**
+ * @file data_type.h
+ * @brief Data storage structures for PRTree
+ *
+ * Contains DataType class for storing index-bounding box pairs
+ * and related utility functions.
+ */
+#pragma once
+
+#include <utility>
+
+#include "prtree/core/detail/bounding_box.h"
+#include "prtree/core/detail/types.h"
+
+// Phase 8: Apply C++20 concept constraints
+template <IndexType T, int D = 2> class DataType {
+public:
+  BB<D> second;
+  T first;
+
+  DataType() noexcept = default;
+
+  DataType(const T &f, const BB<D> &s) {
+    first = f;
+    second = s;
+  }
+
+  DataType(T &&f, BB<D> &&s) noexcept {
+    first = std::move(f);
+    second = std::move(s);
+  }
+
+  void swap(DataType& other) noexcept {
+    using std::swap;
+    swap(first, other.first);
+    swap(second, other.second);
+  }
+
+  template <class Archive> void serialize(Archive &ar) { ar(first, second); }
+};
+
+template <class T, int D = 2>
+void clean_data(DataType<T, D> *b, DataType<T, D> *e) {
+  for (DataType<T, D> *it = e - 1; it >= b; --it) {
+    it->~DataType<T, D>();
+  }
+}