Skip to content

Commit 6548882

Browse files
claudeclaude
committed
Add adaptive epsilon and configurable precision parameters
Addresses Float32 Precision Issues from precision validation gaps: - Add adaptive epsilon calculation that scales with coordinate magnitude - Add configurable precision parameters (relative_epsilon, absolute_epsilon) - Implement subnormal number detection in validate_box() - Update both insert() overloads to use adaptive epsilon - Add precision control methods (setters/getters) to PRTree class - Expose precision control to Python via pybind11 bindings - Add comprehensive test suite for adaptive epsilon behavior This improves precision handling across different coordinate scales, from small (< 1.0) to large (> 1e6) magnitudes. Note: Current architecture still forces float32 tree structure on all users. Float64 data only used for idx2exact refinement. Future work should consider separating float32/float64 builds or templating Real type.
1 parent def408d commit 6548882

File tree

3 files changed

+455
-2
lines changed

3 files changed

+455
-2
lines changed

‎include/prtree/core/prtree.h‎

Lines changed: 121 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,50 @@ template <IndexType T, int B = 6, int D = 2> class PRTree {
7777

7878
mutable std::unique_ptr<std::recursive_mutex> tree_mutex_;
7979

80+
// Precision control parameters
81+
Real relative_epsilon_ = 1e-6; // Relative epsilon for adaptive precision
82+
Real absolute_epsilon_ = 1e-8; // Absolute epsilon (backward compatible default)
83+
bool use_adaptive_epsilon_ = true; // Use adaptive epsilon based on coordinate magnitude
84+
bool detect_subnormal_ = true; // Detect and handle subnormal numbers
85+
86+
// Helper: Calculate adaptive epsilon for insert operations
87+
Real calculate_adaptive_epsilon(const BB<D> &bb) const {
88+
if (!use_adaptive_epsilon_) {
89+
return absolute_epsilon_;
90+
}
91+
92+
// Calculate the maximum extent of the bounding box
93+
Real max_extent = 0.0;
94+
for (int i = 0; i < D; ++i) {
95+
Real extent = bb.max(i) - bb.min(i);
96+
if (extent > max_extent) {
97+
max_extent = extent;
98+
}
99+
}
100+
101+
// For degenerate boxes (points), use the magnitude of coordinates
102+
if (max_extent < std::numeric_limits<Real>::epsilon()) {
103+
Real max_magnitude = 0.0;
104+
for (int i = 0; i < D; ++i) {
105+
Real mag = std::max(std::abs(bb.min(i)), std::abs(bb.max(i)));
106+
if (mag > max_magnitude) {
107+
max_magnitude = mag;
108+
}
109+
}
110+
max_extent = max_magnitude;
111+
}
112+
113+
// Adaptive epsilon: relative to the scale + absolute minimum
114+
// This ensures reasonable behavior across different coordinate scales
115+
Real adaptive_eps = max_extent * relative_epsilon_ + absolute_epsilon_;
116+
117+
// Clamp to reasonable bounds to avoid numerical issues
118+
Real min_epsilon = std::numeric_limits<Real>::epsilon() * 10.0;
119+
Real max_epsilon = max_extent * 0.01; // At most 1% of the extent
120+
121+
return std::max(min_epsilon, std::min(adaptive_eps, max_epsilon));
122+
}
123+
80124
public:
81125
template <class Archive> void serialize(Archive &archive) {
82126
archive(flat_tree, idx2bb, idx2data, global_idx, n_at_build, idx2exact);
@@ -126,6 +170,20 @@ template <IndexType T, int B = 6, int D = 2> class PRTree {
126170
"Bounding box coordinates must be finite (no NaN or Inf)");
127171
}
128172

173+
// Check for subnormal numbers if detection is enabled
174+
if (detect_subnormal_) {
175+
// A number is subnormal if it's non-zero but not normal
176+
bool min_subnormal = (min_val != 0.0) && !std::isnormal(min_val);
177+
bool max_subnormal = (max_val != 0.0) && !std::isnormal(max_val);
178+
179+
if (min_subnormal || max_subnormal) {
180+
throw std::runtime_error(
181+
"Bounding box contains subnormal numbers which may cause "
182+
"precision issues. Consider rescaling coordinates or using "
183+
"larger values. Subnormal detection can be disabled if needed.");
184+
}
185+
}
186+
129187
// Enforce min <= max
130188
if (min_val > max_val) {
131189
throw std::runtime_error(
@@ -357,9 +415,11 @@ template <IndexType T, int B = 6, int D = 2> class PRTree {
357415
idx2bb.emplace(idx, bb);
358416
set_obj(idx, objdumps);
359417

418+
// Use adaptive epsilon based on bounding box scale
419+
Real adaptive_eps = calculate_adaptive_epsilon(bb);
360420
Real delta[D];
361421
for (int i = 0; i < D; ++i) {
362-
delta[i] = bb.max(i) - bb.min(i) + 0.00000001;
422+
delta[i] = bb.max(i) - bb.min(i) + adaptive_eps;
363423
}
364424

365425
// find the leaf node to insert
@@ -503,9 +563,11 @@ template <IndexType T, int B = 6, int D = 2> class PRTree {
503563
idx2exact[idx] = exact_coords; // Store exact coordinates for refinement
504564
set_obj(idx, objdumps);
505565

566+
// Use adaptive epsilon based on bounding box scale
567+
Real adaptive_eps = calculate_adaptive_epsilon(bb);
506568
Real delta[D];
507569
for (int i = 0; i < D; ++i) {
508-
delta[i] = bb.max(i) - bb.min(i) + 0.00000001;
570+
delta[i] = bb.max(i) - bb.min(i) + adaptive_eps;
509571
}
510572

511573
// find the leaf node to insert
@@ -1194,4 +1256,61 @@ template <IndexType T, int B = 6, int D = 2> class PRTree {
11941256
capsule // capsule for cleanup
11951257
);
11961258
}
1259+
1260+
// Precision control methods
1261+
1262+
/**
1263+
* Set relative epsilon for adaptive precision calculation.
1264+
* This epsilon is multiplied by the coordinate scale to determine
1265+
* the precision threshold for insert operations.
1266+
*
1267+
* @param epsilon Relative epsilon (default: 1e-6)
1268+
*/
1269+
void set_relative_epsilon(Real epsilon) {
1270+
if (epsilon <= 0.0 || !std::isfinite(epsilon)) {
1271+
throw std::runtime_error("Relative epsilon must be positive and finite");
1272+
}
1273+
relative_epsilon_ = epsilon;
1274+
}
1275+
1276+
/**
1277+
* Set absolute epsilon for precision calculation.
1278+
* This is the minimum epsilon used regardless of coordinate scale,
1279+
* ensuring backward compatibility and reasonable behavior for small coordinates.
1280+
*
1281+
* @param epsilon Absolute epsilon (default: 1e-8)
1282+
*/
1283+
void set_absolute_epsilon(Real epsilon) {
1284+
if (epsilon <= 0.0 || !std::isfinite(epsilon)) {
1285+
throw std::runtime_error("Absolute epsilon must be positive and finite");
1286+
}
1287+
absolute_epsilon_ = epsilon;
1288+
}
1289+
1290+
/**
1291+
* Enable or disable adaptive epsilon calculation.
1292+
* When disabled, only absolute_epsilon is used (backward compatible behavior).
1293+
*
1294+
* @param enabled True to enable adaptive epsilon (default: true)
1295+
*/
1296+
void set_adaptive_epsilon(bool enabled) {
1297+
use_adaptive_epsilon_ = enabled;
1298+
}
1299+
1300+
/**
1301+
* Enable or disable subnormal number detection.
1302+
* When enabled, validation will reject subnormal numbers to avoid precision issues.
1303+
*
1304+
* @param enabled True to enable detection (default: true)
1305+
*/
1306+
void set_subnormal_detection(bool enabled) {
1307+
detect_subnormal_ = enabled;
1308+
}
1309+
1310+
// Getters for precision parameters
1311+
1312+
Real get_relative_epsilon() const { return relative_epsilon_; }
1313+
Real get_absolute_epsilon() const { return absolute_epsilon_; }
1314+
bool get_adaptive_epsilon() const { return use_adaptive_epsilon_; }
1315+
bool get_subnormal_detection() const { return detect_subnormal_; }
11971316
};

‎src/cpp/bindings/python_bindings.cc‎

Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,42 @@ PYBIND11_MODULE(PRTree, m) {
8080
Find all pairs of intersecting AABBs.
8181
Returns a numpy array of shape (n_pairs, 2) where each row contains
8282
a pair of indices (i, j) with i < j representing intersecting AABBs.
83+
)pbdoc")
84+
.def("set_relative_epsilon", &PRTree<T, B, 2>::set_relative_epsilon,
85+
py::arg("epsilon"),
86+
R"pbdoc(
87+
Set relative epsilon for adaptive precision calculation.
88+
)pbdoc")
89+
.def("set_absolute_epsilon", &PRTree<T, B, 2>::set_absolute_epsilon,
90+
py::arg("epsilon"),
91+
R"pbdoc(
92+
Set absolute epsilon for precision calculation.
93+
)pbdoc")
94+
.def("set_adaptive_epsilon", &PRTree<T, B, 2>::set_adaptive_epsilon,
95+
py::arg("enabled"),
96+
R"pbdoc(
97+
Enable or disable adaptive epsilon calculation.
98+
)pbdoc")
99+
.def("set_subnormal_detection", &PRTree<T, B, 2>::set_subnormal_detection,
100+
py::arg("enabled"),
101+
R"pbdoc(
102+
Enable or disable subnormal number detection.
103+
)pbdoc")
104+
.def("get_relative_epsilon", &PRTree<T, B, 2>::get_relative_epsilon,
105+
R"pbdoc(
106+
Get current relative epsilon value.
107+
)pbdoc")
108+
.def("get_absolute_epsilon", &PRTree<T, B, 2>::get_absolute_epsilon,
109+
R"pbdoc(
110+
Get current absolute epsilon value.
111+
)pbdoc")
112+
.def("get_adaptive_epsilon", &PRTree<T, B, 2>::get_adaptive_epsilon,
113+
R"pbdoc(
114+
Check if adaptive epsilon is enabled.
115+
)pbdoc")
116+
.def("get_subnormal_detection", &PRTree<T, B, 2>::get_subnormal_detection,
117+
R"pbdoc(
118+
Check if subnormal detection is enabled.
83119
)pbdoc");
84120

85121
py::class_<PRTree<T, B, 3>>(m, "_PRTree3D")
@@ -146,6 +182,42 @@ PYBIND11_MODULE(PRTree, m) {
146182
Find all pairs of intersecting AABBs.
147183
Returns a numpy array of shape (n_pairs, 2) where each row contains
148184
a pair of indices (i, j) with i < j representing intersecting AABBs.
185+
)pbdoc")
186+
.def("set_relative_epsilon", &PRTree<T, B, 3>::set_relative_epsilon,
187+
py::arg("epsilon"),
188+
R"pbdoc(
189+
Set relative epsilon for adaptive precision calculation.
190+
)pbdoc")
191+
.def("set_absolute_epsilon", &PRTree<T, B, 3>::set_absolute_epsilon,
192+
py::arg("epsilon"),
193+
R"pbdoc(
194+
Set absolute epsilon for precision calculation.
195+
)pbdoc")
196+
.def("set_adaptive_epsilon", &PRTree<T, B, 3>::set_adaptive_epsilon,
197+
py::arg("enabled"),
198+
R"pbdoc(
199+
Enable or disable adaptive epsilon calculation.
200+
)pbdoc")
201+
.def("set_subnormal_detection", &PRTree<T, B, 3>::set_subnormal_detection,
202+
py::arg("enabled"),
203+
R"pbdoc(
204+
Enable or disable subnormal number detection.
205+
)pbdoc")
206+
.def("get_relative_epsilon", &PRTree<T, B, 3>::get_relative_epsilon,
207+
R"pbdoc(
208+
Get current relative epsilon value.
209+
)pbdoc")
210+
.def("get_absolute_epsilon", &PRTree<T, B, 3>::get_absolute_epsilon,
211+
R"pbdoc(
212+
Get current absolute epsilon value.
213+
)pbdoc")
214+
.def("get_adaptive_epsilon", &PRTree<T, B, 3>::get_adaptive_epsilon,
215+
R"pbdoc(
216+
Check if adaptive epsilon is enabled.
217+
)pbdoc")
218+
.def("get_subnormal_detection", &PRTree<T, B, 3>::get_subnormal_detection,
219+
R"pbdoc(
220+
Check if subnormal detection is enabled.
149221
)pbdoc");
150222

151223
py::class_<PRTree<T, B, 4>>(m, "_PRTree4D")
@@ -212,6 +284,42 @@ PYBIND11_MODULE(PRTree, m) {
212284
Find all pairs of intersecting AABBs.
213285
Returns a numpy array of shape (n_pairs, 2) where each row contains
214286
a pair of indices (i, j) with i < j representing intersecting AABBs.
287+
)pbdoc")
288+
.def("set_relative_epsilon", &PRTree<T, B, 4>::set_relative_epsilon,
289+
py::arg("epsilon"),
290+
R"pbdoc(
291+
Set relative epsilon for adaptive precision calculation.
292+
)pbdoc")
293+
.def("set_absolute_epsilon", &PRTree<T, B, 4>::set_absolute_epsilon,
294+
py::arg("epsilon"),
295+
R"pbdoc(
296+
Set absolute epsilon for precision calculation.
297+
)pbdoc")
298+
.def("set_adaptive_epsilon", &PRTree<T, B, 4>::set_adaptive_epsilon,
299+
py::arg("enabled"),
300+
R"pbdoc(
301+
Enable or disable adaptive epsilon calculation.
302+
)pbdoc")
303+
.def("set_subnormal_detection", &PRTree<T, B, 4>::set_subnormal_detection,
304+
py::arg("enabled"),
305+
R"pbdoc(
306+
Enable or disable subnormal number detection.
307+
)pbdoc")
308+
.def("get_relative_epsilon", &PRTree<T, B, 4>::get_relative_epsilon,
309+
R"pbdoc(
310+
Get current relative epsilon value.
311+
)pbdoc")
312+
.def("get_absolute_epsilon", &PRTree<T, B, 4>::get_absolute_epsilon,
313+
R"pbdoc(
314+
Get current absolute epsilon value.
315+
)pbdoc")
316+
.def("get_adaptive_epsilon", &PRTree<T, B, 4>::get_adaptive_epsilon,
317+
R"pbdoc(
318+
Check if adaptive epsilon is enabled.
319+
)pbdoc")
320+
.def("get_subnormal_detection", &PRTree<T, B, 4>::get_subnormal_detection,
321+
R"pbdoc(
322+
Check if subnormal detection is enabled.
215323
)pbdoc");
216324

217325
#ifdef VERSION_INFO

0 commit comments

Comments
 (0)