From a380b83518f3ecd195af9bf9eae324840777d5db Mon Sep 17 00:00:00 2001 From: Riddho Ridwanul Haque <33618344+RiddhoHaque@users.noreply.github.com> Date: Thu, 7 May 2026 15:58:12 -0400 Subject: [PATCH 1/3] Add incremental miniball updates Add C++ warm-start constructor and contains helper for incremental point insertion. Expose incremental_miniball in the Python API. Add regression tests comparing incremental updates against full recomputation. Fix Windows build compatibility in Seb_debug by avoiding Unix-only timer headers on MSVC. --- cpp/main/Seb-inl.h | 152 ++++++++++++++++++++++++++++++++--- cpp/main/Seb.h | 43 +++++++++- cpp/main/Seb_debug.C | 16 ++++ cpp/main/Seb_debug.h | 8 ++ cpp/main/Subspan-inl.h | 19 +++++ cpp/main/Subspan.h | 2 + python/miniball/__init__.py | 66 ++++++++++++++- python/miniball_python.cpp | 120 +++++++++++++++++++++++++-- python/test/test_miniball.py | 87 +++++++++++++++++++- 9 files changed, 489 insertions(+), 24 deletions(-) diff --git a/cpp/main/Seb-inl.h b/cpp/main/Seb-inl.h index 3aa920b..3c64e53 100644 --- a/cpp/main/Seb-inl.h +++ b/cpp/main/Seb-inl.h @@ -106,6 +106,43 @@ namespace SEB_NAMESPACE { return false; } + template + Float Smallest_enclosing_ball::find_stop_fraction( + const Float* direction, int& stopper) + { + using std::inner_product; + + Float scale = 0; + stopper = -1; + const Pt& support_point = S[support->any_member()]; + + for (unsigned int j = 0; j < S.size(); ++j) + if (!support->is_member(j)) { + Float dist = 0; + for (unsigned int i = 0; i < dim; ++i) { + center_to_point[i] = S[j][i] - center[i]; + dist += sqr(center_to_point[i]); + } + + for (unsigned int i = 0; i < dim; ++i) + center_to_point[i] = S[j][i] - support_point[i]; + const Float denom = + 2 * inner_product(direction, direction+dim, + center_to_point, Float(0)); + if (denom == 0) + continue; + + const Float bound = (dist - radius_square) / denom; + if (bound > 0 && (stopper < 0 || bound < scale)) { + scale = bound; + stopper = j; + } + } + + return scale; + } + + template Float Smallest_enclosing_ball::find_stop_fraction(int& stopper) // Given the center of the current enclosing ball and the @@ -173,7 +210,7 @@ namespace SEB_NAMESPACE { template - void Smallest_enclosing_ball::update() + void Smallest_enclosing_ball::pivot() // The main function containing the main loop. // Iteratively, we compute the point in support that is closest // to the current center and then walk towards this target as far @@ -187,17 +224,6 @@ namespace SEB_NAMESPACE { { SEB_DEBUG (int iteration = 0;) - SEB_TIMER_START("computation"); - - // optimistically, we set this flag now; - // on return from this function it will be true: - up_to_date = true; - - init_ball(); - - // Invariant: The ball B(center,radius_) always contains the whole - // point set S and has the points in support on its boundary. - while (true) { SEB_LOG ("debug"," iteration " << ++iteration << std::endl); @@ -293,6 +319,108 @@ namespace SEB_NAMESPACE { } } + template + void Smallest_enclosing_ball::update( + const Float* previous_center, Float previous_squared_radius, + unsigned int new_point_index) + { + SEB_TIMER_START("computation"); + + // optimistically, we set this flag now; + // on return from this function it will be true: + up_to_date = true; + + // Incremental construction starts from the smallest ball containing the + // previous miniball and the new outside point. + if (previous_center != NULL) { + SEB_ASSERT(S.size() > 0); + SEB_ASSERT(new_point_index < S.size()); + const Pt& new_point = S[new_point_index]; + Float dist_square = 0; + for (unsigned int i = 0; i < dim; ++i) + dist_square += sqr(new_point[i] - previous_center[i]); + const Float previous_radius = sqrt(previous_squared_radius); + const Float dist = sqrt(dist_square); + const Float new_radius = (dist + previous_radius) / 2; + const Float shift = (dist - previous_radius) / (2 * dist); + for (unsigned int i = 0; i < dim; ++i) + center[i] = previous_center[i] + shift * (new_point[i] - previous_center[i]); + radius_square = sqr(new_radius); + radius_ = new_radius; + if (support != NULL) + support->reset(new_point_index); + else + support = new Subspan(dim, S, new_point_index); + SEB_STATS(entry_count = std::vector(S.size(),0)); + } else { + init_ball(); + } + + // Invariant: The ball B(center,radius_) always contains the whole + // point set S and has the points in support on its boundary. + + pivot(); + } + + template + void Smallest_enclosing_ball::append_point( + unsigned int new_point_index) + { + if (!up_to_date) + update(); + + SEB_ASSERT(new_point_index < S.size()); + + if (support != NULL) + support->resize_membership(); + + Float dist = 0; + for (unsigned int i = 0; i < dim; ++i) + dist += sqr(S[new_point_index][i] - center[i]); + if (dist <= radius_square) + return; + + while (!contains(S[new_point_index])) { + if (support->size() > dim) { + update(); + return; + } + + support->shortest_vector_to_span(center, center_to_aff); + support->shortest_vector_to_span(S[new_point_index], center_to_point); + dist_to_aff_square = 0; + for (unsigned int i = 0; i < dim; ++i) { + center_to_aff[i] -= center_to_point[i]; + dist_to_aff_square += sqr(center_to_aff[i]); + } + dist_to_aff = sqrt(dist_to_aff_square); + if (dist_to_aff <= Eps * radius_) { + update(); + return; + } + + int stopper; + Float scale = find_stop_fraction(center_to_aff, stopper); + if (stopper < 0) { + update(); + return; + } + + for (unsigned int i = 0; i < dim; ++i) + center[i] += scale * center_to_aff[i]; + + const Pt& stop_point = S[support->any_member()]; + radius_square = 0; + for (unsigned int i = 0; i < dim; ++i) + radius_square += sqr(stop_point[i] - center[i]); + radius_ = sqrt(radius_square); + + support->add_point(stopper); + } + + pivot(); + } + template void Smallest_enclosing_ball::verify() { diff --git a/cpp/main/Seb.h b/cpp/main/Seb.h index 895de03..b49c5ad 100644 --- a/cpp/main/Seb.h +++ b/cpp/main/Seb.h @@ -42,6 +42,27 @@ namespace SEB_NAMESPACE { update(); } + Smallest_enclosing_ball(unsigned int d, const PointAccessor &P, + const Float* previous_center, + Float previous_squared_radius, + unsigned int new_point_index) + // Constructs an instance representing the miniball of points from + // set S, using the miniball of S without new_point_index as a warm + // start. The new point is assumed to lie outside that previous ball. + : dim(d), S(P), up_to_date(true), support(NULL) + { + allocate_resources(); + SEB_ASSERT(!is_empty()); + update(previous_center, previous_squared_radius, new_point_index); + // Warm starting is an optimization; keep the constructor exact if that + // start does not preserve the enclosing-ball invariant. + for (unsigned int j = 0; j < S.size(); ++j) + if (!contains(S[j])) { + update(); + break; + } + } + ~Smallest_enclosing_ball() { deallocate_resources(); @@ -58,6 +79,8 @@ namespace SEB_NAMESPACE { up_to_date = false; } + void append_point(unsigned int new_point_index); + public: // access: bool is_empty() @@ -114,6 +137,20 @@ namespace SEB_NAMESPACE { return center+dim; } + bool contains(const Pt& point) + // Returns whether point is contained in the miniball. + // Precondition: !is_empty() + { + if (!up_to_date) + update(); + + SEB_ASSERT(!is_empty()); + Float dist = 0; + for (unsigned int i = 0; i < dim; ++i) + dist += sqr(point[i] - center[i]); + return dist <= radius_square; + } + public: // testing: void verify(); @@ -137,9 +174,13 @@ namespace SEB_NAMESPACE { private: // internal helper routines for the actual algorithm: void init_ball(); Float find_stop_fraction(int& hinderer); + Float find_stop_fraction(const Float* direction, int& hinderer); bool successful_drop(); + void pivot(); - void update(); + void update(const Float* previous_center = NULL, + Float previous_squared_radius = 0, + unsigned int new_point_index = 0); private: // we forbid copying (since we have dynamic storage): Smallest_enclosing_ball(const Smallest_enclosing_ball&); diff --git a/cpp/main/Seb_debug.C b/cpp/main/Seb_debug.C index 6924edf..1199ab0 100644 --- a/cpp/main/Seb_debug.C +++ b/cpp/main/Seb_debug.C @@ -4,8 +4,12 @@ // Kaspar Fischer #include +#ifdef _WIN32 +#include +#else #include #include +#endif #include "Seb_configure.h" @@ -70,6 +74,7 @@ namespace SEB_NAMESPACE { // long tv_usec; /* microseconds */ // }; // +#ifndef _WIN32 inline timeval& operator-=(timeval &t1,const timeval &t2) { t1.tv_sec -= t2.tv_sec; @@ -79,6 +84,7 @@ namespace SEB_NAMESPACE { } return t1; } +#endif Timer::Timer() { @@ -96,6 +102,9 @@ namespace SEB_NAMESPACE { void Timer::start(const char *timer_name) { +#ifdef _WIN32 + timers[std::string(timer_name)] = std::chrono::steady_clock::now(); +#else // fetch current usage: rusage now; int status = getrusage(RUSAGE_SELF,&now); @@ -103,6 +112,7 @@ namespace SEB_NAMESPACE { // save it: timers[std::string(timer_name)] = now.ru_utime; +#endif } float Timer::lapse(const char *name) @@ -110,6 +120,11 @@ namespace SEB_NAMESPACE { // assert that start(name) has been called before: SEB_ASSERT(timers.find(std::string(name)) != timers.end()); +#ifdef _WIN32 + const std::chrono::duration elapsed = + std::chrono::steady_clock::now() - (*timers.find(std::string(name))).second; + return elapsed.count(); +#else // get current usage: rusage now; int status = getrusage(RUSAGE_SELF,&now); @@ -118,6 +133,7 @@ namespace SEB_NAMESPACE { // compute elapsed usage: now.ru_utime -= (*timers.find(std::string(name))).second; return now.ru_utime.tv_sec + now.ru_utime.tv_usec * 1e-6; +#endif } } // namespace SEB_NAMESPACE diff --git a/cpp/main/Seb_debug.h b/cpp/main/Seb_debug.h index c54bd24..ac1bd5e 100644 --- a/cpp/main/Seb_debug.h +++ b/cpp/main/Seb_debug.h @@ -9,8 +9,12 @@ #include #include #include +#ifdef _WIN32 +#include +#else #include #include +#endif namespace SEB_NAMESPACE { @@ -77,7 +81,11 @@ namespace SEB_NAMESPACE { float lapse(const char *name); private: // private members: +#ifdef _WIN32 + typedef std::map Timers; +#else typedef std::map Timers; +#endif Timers timers; // a collection of pairs (k,v) where // k is the timer name and v is the // (started) timer associated with k diff --git a/cpp/main/Subspan-inl.h b/cpp/main/Subspan-inl.h index 404b48f..d76fd89 100644 --- a/cpp/main/Subspan-inl.h +++ b/cpp/main/Subspan-inl.h @@ -88,6 +88,25 @@ namespace SEB_NAMESPACE { delete[] w; } + template + void Subspan::resize_membership() + { + membership.resize(S.size(), false); + } + + template + void Subspan::reset(unsigned int index) + { + membership.assign(S.size(), false); + for (unsigned int i=0; i void Subspan::add_point(int index) { SEB_ASSERT(!is_member(index)); diff --git a/cpp/main/Subspan.h b/cpp/main/Subspan.h index b45cb90..5690161 100644 --- a/cpp/main/Subspan.h +++ b/cpp/main/Subspan.h @@ -89,6 +89,8 @@ namespace SEB_NAMESPACE { void add_point(int global_index); void remove_point(unsigned int local_index); + void resize_membership(); + void reset(unsigned int global_index); public: // access: diff --git a/python/miniball/__init__.py b/python/miniball/__init__.py index 787679b..2af671a 100644 --- a/python/miniball/__init__.py +++ b/python/miniball/__init__.py @@ -1,8 +1,12 @@ -from ._miniball import _compute_miniball +from ._miniball import ( + _Miniball, + _compute_miniball, + _compute_miniball_incremental, +) import numpy as np -__all__ = ["miniball"] +__all__ = ["miniball", "incremental_miniball", "Miniball"] def miniball(points: np.typing.ArrayLike): @@ -14,3 +18,61 @@ def miniball(points: np.typing.ArrayLike): msg = f"Input array must be 2-dimensional! Got shape `{points.shape}`." raise TypeError(msg) return _compute_miniball(points) + + +def incremental_miniball( + points: np.typing.ArrayLike, + current: dict, + point: np.typing.ArrayLike, +): + """Compute the miniball after appending one point to an existing point set.""" + + points = np.ascontiguousarray(points, dtype=np.float64) + if len(points.shape) != 2: + msg = f"Input array must be 2-dimensional! Got shape `{points.shape}`." + raise TypeError(msg) + + point = np.ascontiguousarray(point, dtype=np.float64) + if len(point.shape) != 1: + msg = f"New point must be 1-dimensional! Got shape `{point.shape}`." + raise TypeError(msg) + + center = np.ascontiguousarray(current["center"], dtype=np.float64) + radius_squared = float(current["radius_squared"]) + delta = point - center + distance_squared = float(np.dot(delta, delta)) + + # If the new point is already enclosed, the previous miniball is unchanged. + if distance_squared <= radius_squared: + return current + + # Otherwise the native helper starts from a ball with point on the boundary. + return _compute_miniball_incremental(points, center, radius_squared, point) + + +class Miniball: + """Maintain a miniball as points are appended.""" + + def __init__(self, points: np.typing.ArrayLike): + points = np.ascontiguousarray(points, dtype=np.float64) + if len(points.shape) != 2: + msg = f"Input array must be 2-dimensional! Got shape `{points.shape}`." + raise TypeError(msg) + self._miniball = _Miniball(points) + + def add(self, point: np.typing.ArrayLike): + point = np.ascontiguousarray(point, dtype=np.float64) + if len(point.shape) != 1: + msg = f"New point must be 1-dimensional! Got shape `{point.shape}`." + raise TypeError(msg) + return self._miniball.add(point) + + def add_points(self, points: np.typing.ArrayLike): + points = np.ascontiguousarray(points, dtype=np.float64) + if len(points.shape) != 2: + msg = f"Input array must be 2-dimensional! Got shape `{points.shape}`." + raise TypeError(msg) + return self._miniball.add_points(points) + + def result(self): + return self._miniball.result() diff --git a/python/miniball_python.cpp b/python/miniball_python.cpp index 39bf4a0..27404a9 100644 --- a/python/miniball_python.cpp +++ b/python/miniball_python.cpp @@ -7,13 +7,77 @@ */ #include "../cpp/main/Seb.h" +#include +#include #include #include namespace nb = nanobind; using Point = Seb::Point; -using Miniball = Seb::Smallest_enclosing_ball; +using NativeMiniball = Seb::Smallest_enclosing_ball; + +nb::dict miniball_result(NativeMiniball &mb, size_t dim) { + nb::dict result; + result["center"] = + nb::ndarray(mb.center_begin(), {dim}).cast(); + result["radius"] = mb.radius(); + result["radius_squared"] = mb.squared_radius(); + + return result; +} + +class Miniball { +public: + Miniball( + nb::ndarray, nb::c_contig> points_arr) { + size_t n_points = points_arr.shape(0); + dim_ = points_arr.shape(1); + if (n_points == 0) { + throw std::invalid_argument("Input array must contain at least one point."); + } + + const double *data = points_arr.data(); + points_.reset(new std::vector()); + points_->reserve(n_points); + for (size_t i = 0; i < n_points; ++i) { + points_->emplace_back(dim_, data + i * dim_); + } + miniball_.reset(new NativeMiniball(static_cast(dim_), *points_)); + } + + nb::dict add(nb::ndarray, nb::c_contig> point_arr) { + if (point_arr.shape(0) != dim_) { + throw std::invalid_argument("Point must match the point-set dimension."); + } + add_point(point_arr.data()); + return result(); + } + + nb::dict add_points( + nb::ndarray, nb::c_contig> points_arr) { + if (points_arr.shape(1) != dim_) { + throw std::invalid_argument("Points must match the point-set dimension."); + } + const double *data = points_arr.data(); + for (size_t i = 0; i < points_arr.shape(0); ++i) { + add_point(data + i * dim_); + } + return result(); + } + + nb::dict result() { return miniball_result(*miniball_, dim_); } + +private: + void add_point(const double *point) { + points_->emplace_back(dim_, point); + miniball_->append_point(static_cast(points_->size() - 1)); + } + + size_t dim_; + std::unique_ptr> points_; + std::unique_ptr miniball_; +}; /** * @brief Computes the smallest enclosing ball for a set of points. @@ -42,15 +106,46 @@ nb::dict compute_miniball( } // Compute the smallest enclosing ball. - Miniball mb(dim, points); + NativeMiniball mb(dim, points); - nb::dict result; - result["center"] = - nb::ndarray(mb.center_begin(), {dim}).cast(); - result["radius"] = mb.radius(); - result["radius_squared"] = mb.squared_radius(); + return miniball_result(mb, dim); +} - return result; +/** + * @brief Computes the smallest enclosing ball after appending one point. + * + * The caller supplies the existing point set, its current miniball, and a new + * point outside that ball. The search starts from the smallest ball that + * contains the old ball and has the new point on its boundary. + */ +nb::dict compute_miniball_incremental( + nb::ndarray, nb::c_contig> points_arr, + nb::ndarray, nb::c_contig> center_arr, + double radius_squared, + nb::ndarray, nb::c_contig> point_arr) { + size_t n_points = points_arr.shape(0); + size_t dim = points_arr.shape(1); + + if (center_arr.shape(0) != dim || point_arr.shape(0) != dim) { + throw std::invalid_argument( + "Center and point must match the point-set dimension."); + } + + const double *data = points_arr.data(); + const double *old_center = center_arr.data(); + const double *new_point = point_arr.data(); + + std::vector points; + points.reserve(n_points + 1); + for (size_t i = 0; i < n_points; ++i) { + points.emplace_back(dim, data + i * dim); + } + points.emplace_back(dim, new_point); + + NativeMiniball mb(dim, points, old_center, radius_squared, + static_cast(n_points)); + + return miniball_result(mb, dim); } // Define the Python module using the NB_MODULE macro. @@ -58,4 +153,13 @@ nb::dict compute_miniball( NB_MODULE(_miniball, m) { m.def("_compute_miniball", &compute_miniball, nb::arg("points"), "Compute the smallest enclosing ball for a set of points."); + m.def("_compute_miniball_incremental", &compute_miniball_incremental, + nb::arg("points"), nb::arg("center"), nb::arg("radius_squared"), + nb::arg("point"), + "Compute the smallest enclosing ball after appending one outside point."); + nb::class_(m, "_Miniball") + .def(nb::init, nb::c_contig>>()) + .def("add", &Miniball::add) + .def("add_points", &Miniball::add_points) + .def("result", &Miniball::result); } diff --git a/python/test/test_miniball.py b/python/test/test_miniball.py index ac42bc9..b59ae6e 100644 --- a/python/test/test_miniball.py +++ b/python/test/test_miniball.py @@ -1,5 +1,5 @@ import pytest -from miniball import miniball +from miniball import Miniball, incremental_miniball, miniball import numpy as np @@ -37,3 +37,88 @@ def test_three_points(): res = miniball(test_vector) np.testing.assert_allclose(res["center"], [0.5, 0.5]) assert res["radius_squared"] == 0.5 + + +def test_incremental_point_inside_current_ball(): + test_vector = np.array([[0, 0], [2, 0]], dtype=np.double) + current = miniball(test_vector) + + res = incremental_miniball(test_vector, current, [1, 0]) + + assert res is current + + +def test_incremental_point_outside_matches_full_recomputation(): + test_vector = np.array([[0, 0], [2, 0], [1, 1]], dtype=np.double) + point = np.array([1, 4], dtype=np.double) + current = miniball(test_vector) + + incremental = incremental_miniball(test_vector, current, point) + full = miniball(np.vstack([test_vector, point])) + + np.testing.assert_allclose(incremental["center"], full["center"]) + assert incremental["radius"] == pytest.approx(full["radius"]) + assert incremental["radius_squared"] == pytest.approx(full["radius_squared"]) + + +def test_incremental_random_high_dimensional_sequence_matches_full_recomputation(): + rng = np.random.default_rng(0) + points = rng.normal(size=(2000, 100)) + current = miniball(points) + + for point in rng.normal(size=(100, 100)): + previous = points + points = np.vstack([points, point]) + current = incremental_miniball(previous, current, point) + full = miniball(points) + + np.testing.assert_allclose(current["center"], full["center"], rtol=1e-10) + assert current["radius"] == pytest.approx(full["radius"], rel=1e-10) + assert current["radius_squared"] == pytest.approx( + full["radius_squared"], rel=1e-10 + ) + + +def test_stateful_incremental_result_matches_full_recomputation(): + points = np.array([[0, 0], [2, 0], [1, 1]], dtype=np.double) + stream = Miniball(points) + + np.testing.assert_allclose(stream.result()["center"], miniball(points)["center"]) + + for point in np.array([[1, 4], [-1, 2], [3, 3]], dtype=np.double): + points = np.vstack([points, point]) + current = stream.add(point) + full = miniball(points) + + np.testing.assert_allclose(current["center"], full["center"]) + assert current["radius"] == pytest.approx(full["radius"]) + assert current["radius_squared"] == pytest.approx(full["radius_squared"]) + + +def test_stateful_incremental_add_points_matches_repeated_add(): + initial = np.array([[0, 0], [2, 0], [1, 1]], dtype=np.double) + additions = np.array([[1, 4], [-1, 2], [3, 3]], dtype=np.double) + repeated = Miniball(initial) + batched = Miniball(initial) + + for point in additions: + repeated_result = repeated.add(point) + batched_result = batched.add_points(additions) + + np.testing.assert_allclose(batched_result["center"], repeated_result["center"]) + assert batched_result["radius"] == pytest.approx(repeated_result["radius"]) + assert batched_result["radius_squared"] == pytest.approx( + repeated_result["radius_squared"] + ) + + +def test_stateful_incremental_dimension_mismatch(): + stream = Miniball(np.array([[0, 0], [2, 0]], dtype=np.double)) + + with pytest.raises(Exception, match="dimension"): + stream.add([1, 2, 3]) + + +def test_stateful_incremental_rejects_empty_initial_points(): + with pytest.raises(Exception, match="at least one point"): + Miniball(np.empty((0, 2), dtype=np.double)) From e3364f74da40ad61ffcdc7ef76508c83e4b2d4f0 Mon Sep 17 00:00:00 2001 From: Riddho Ridwanul Haque <33618344+RiddhoHaque@users.noreply.github.com> Date: Mon, 11 May 2026 12:34:09 -0400 Subject: [PATCH 2/3] Benchmarking Incremental Updates to the Miniball --- cpp/main/Seb-inl.h | 10 +- cpp/main/Seb.h | 13 +- python/benchmark_incremental.py | 277 ++++++++++++++++++++++++++++++++ python/miniball_python.cpp | 1 + 4 files changed, 296 insertions(+), 5 deletions(-) create mode 100644 python/benchmark_incremental.py diff --git a/cpp/main/Seb-inl.h b/cpp/main/Seb-inl.h index 3c64e53..98c74dd 100644 --- a/cpp/main/Seb-inl.h +++ b/cpp/main/Seb-inl.h @@ -222,11 +222,10 @@ namespace SEB_NAMESPACE { // If such an attempt to drop fails, we are done; because then // the center lies even conv(support). { - SEB_DEBUG (int iteration = 0;) - while (true) { - SEB_LOG ("debug"," iteration " << ++iteration << std::endl); + ++last_iteration_count; + SEB_LOG ("debug"," iteration " << last_iteration_count << std::endl); SEB_LOG ("debug"," " << support->size() << " points on boundary" << std::endl); @@ -325,6 +324,7 @@ namespace SEB_NAMESPACE { unsigned int new_point_index) { SEB_TIMER_START("computation"); + last_iteration_count = 0; // optimistically, we set this flag now; // on return from this function it will be true: @@ -369,6 +369,8 @@ namespace SEB_NAMESPACE { if (!up_to_date) update(); + last_iteration_count = 0; + SEB_ASSERT(new_point_index < S.size()); if (support != NULL) @@ -381,6 +383,8 @@ namespace SEB_NAMESPACE { return; while (!contains(S[new_point_index])) { + ++last_iteration_count; + if (support->size() > dim) { update(); return; diff --git a/cpp/main/Seb.h b/cpp/main/Seb.h index b49c5ad..db8241e 100644 --- a/cpp/main/Seb.h +++ b/cpp/main/Seb.h @@ -35,7 +35,8 @@ namespace SEB_NAMESPACE { // Constructs an instance representing the miniball of points from // set S. The dimension of the ambient space is fixed to d for // lifetime of the instance. - : dim(d), S(P), up_to_date(true), support(NULL) + : dim(d), S(P), up_to_date(true), support(NULL), + last_iteration_count(0) { allocate_resources(); SEB_ASSERT(!is_empty()); @@ -49,7 +50,8 @@ namespace SEB_NAMESPACE { // Constructs an instance representing the miniball of points from // set S, using the miniball of S without new_point_index as a warm // start. The new point is assumed to lie outside that previous ball. - : dim(d), S(P), up_to_date(true), support(NULL) + : dim(d), S(P), up_to_date(true), support(NULL), + last_iteration_count(0) { allocate_resources(); SEB_ASSERT(!is_empty()); @@ -151,6 +153,12 @@ namespace SEB_NAMESPACE { return dist <= radius_square; } + unsigned int iterations() + // Returns the number of iterations used by the most recent update. + { + return last_iteration_count; + } + public: // testing: void verify(); @@ -202,6 +210,7 @@ namespace SEB_NAMESPACE { Float *center_to_point; Float *lambdas; Float dist_to_aff, dist_to_aff_square; + unsigned int last_iteration_count; #ifdef SEB_STATS_MODE private: // memeber fields for statistics diff --git a/python/benchmark_incremental.py b/python/benchmark_incremental.py new file mode 100644 index 0000000..db40abd --- /dev/null +++ b/python/benchmark_incremental.py @@ -0,0 +1,277 @@ +from __future__ import annotations + +import argparse +import csv +import statistics +import time +from pathlib import Path + +import numpy as np + +from miniball import Miniball, miniball + + +def unit_direction(rng: np.random.Generator, dim: int) -> np.ndarray: + direction = rng.normal(size=dim) + norm = np.linalg.norm(direction) + if norm == 0: + direction[0] = 1.0 + return direction + return direction / norm + + +def generate_outside_sequence( + rng: np.random.Generator, + initial_points: np.ndarray, + count: int, + distance_factor: float, +) -> np.ndarray: + points = np.asarray(initial_points, dtype=np.float64) + stream = Miniball(points) + current = stream.result() + additions = [] + + for _ in range(count): + center = np.asarray(current["center"], dtype=np.float64) + radius = float(current["radius"]) + distance = distance_factor * radius if radius > 0 else distance_factor + point = center + distance * unit_direction(rng, points.shape[1]) + additions.append(point) + + points = np.vstack([points, point]) + current = stream.add(point) + + return np.asarray(additions, dtype=np.float64) + + +def time_full_recomputation( + all_points: np.ndarray, initial_count: int +) -> tuple[float, list[np.ndarray], dict[str, float]]: + centers = [] + iterations = [] + start = time.perf_counter() + for size in range(initial_count + 1, len(all_points) + 1): + current = miniball(all_points[:size]) + centers.append(np.asarray(current["center"])) + iterations.append(int(current["iterations"])) + elapsed = time.perf_counter() - start + return elapsed, centers, iteration_summary(iterations, "full") + + +def time_incremental( + all_points: np.ndarray, initial_count: int, centers: list[np.ndarray] +) -> tuple[float, dict[str, float]]: + stream = Miniball(all_points[:initial_count]) + iterations = [] + + start = time.perf_counter() + for index in range(initial_count, len(all_points)): + current = stream.add(all_points[index]) + iterations.append(int(current["iterations"])) + np.testing.assert_allclose( + current["center"], + centers[index - initial_count], + rtol=1e-10, + atol=1e-10, + ) + elapsed = time.perf_counter() - start + return elapsed, iteration_summary(iterations, "update") + + +def iteration_summary(iterations: list[int], prefix: str) -> dict[str, float]: + return { + f"{prefix}_max_iterations": max(iterations) if iterations else 0, + f"{prefix}_median_iterations": ( + statistics.median(iterations) if iterations else 0 + ), + f"{prefix}_average_iterations": ( + statistics.fmean(iterations) if iterations else 0.0 + ), + } + + +def run_case( + rng: np.random.Generator, + dim: int, + additions_count: int, + initial_count: int, + distance_factor: float, +) -> dict[str, float]: + initial_points = rng.normal(size=(initial_count, dim)) + additions = generate_outside_sequence( + rng, initial_points, additions_count, distance_factor + ) + all_points = np.vstack([initial_points, additions]) + + full_seconds, centers, full_iteration_stats = time_full_recomputation( + all_points, initial_count + ) + print('Finished full recomputation for dimension {}, point insertions {}'.format(dim, additions_count)) + print('Full recomputation took {:.4f} seconds'.format(full_seconds)) + print( + "Full recomputation iterations: max {}, median {:.1f}, average {:.2f}".format( + full_iteration_stats["full_max_iterations"], + full_iteration_stats["full_median_iterations"], + full_iteration_stats["full_average_iterations"], + ) + ) + update_seconds, iteration_stats = time_incremental( + all_points, initial_count, centers + ) + print('Finished updates for dimension {}, point insertions {}'.format(dim, additions_count)) + print('Updates took {:.4f} seconds'.format(update_seconds)) + print( + "Iterations per update: max {}, median {:.1f}, average {:.2f}".format( + iteration_stats["update_max_iterations"], + iteration_stats["update_median_iterations"], + iteration_stats["update_average_iterations"], + ) + ) + return { + "dimension": dim, + "additions": additions_count, + "full_seconds": full_seconds, + "update_seconds": update_seconds, + "speedup": full_seconds / update_seconds + if update_seconds > 0 + else float("inf"), + **full_iteration_stats, + **iteration_stats, + } + + +def write_csv(path: Path, rows: list[dict[str, float]]) -> None: + with path.open("w", newline="") as file: + writer = csv.DictWriter(file, fieldnames=list(rows[0].keys())) + writer.writeheader() + writer.writerows(rows) + + +def print_table(title: str, rows: list[dict[str, float]], x_key: str) -> None: + print(f"\n{title}") + print("-" * len(title)) + print( + f"{x_key:>10} {'full (s)':>12} {'updates (s)':>17} " + f"{'speedup':>10} {'full max':>8} {'full med':>8} " + f"{'full avg':>8} {'upd max':>8} {'upd med':>8} {'upd avg':>8}" + ) + for row in rows: + print( + f"{int(row[x_key]):>10} " + f"{row['full_seconds']:>12.4f} " + f"{row['update_seconds']:>17.4f} " + f"{row['speedup']:>10.2f} " + f"{row['full_max_iterations']:>8.0f} " + f"{row['full_median_iterations']:>8.1f} " + f"{row['full_average_iterations']:>8.2f} " + f"{row['update_max_iterations']:>8.0f} " + f"{row['update_median_iterations']:>8.1f} " + f"{row['update_average_iterations']:>8.2f}" + ) + + +def plot_results(path: Path, title: str, x_label: str, x_values, rows) -> None: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + full = [row["full_seconds"] for row in rows] + updates = [row["update_seconds"] for row in rows] + + fig, ax = plt.subplots(figsize=(8, 5)) + ax.plot(x_values, full, marker="o", label="Full Recomputation") + ax.plot(x_values, updates, marker="o", label="Warm Started Balls") + ax.set_title(title) + ax.set_xlabel(x_label) + ax.set_ylabel("Runtime (seconds)") + ax.grid(True, which="both", alpha=0.3) + ax.legend() + fig.tight_layout() + fig.savefig(path, dpi=200) + plt.close(fig) + + +def plot_iteration_results( + path: Path, title: str, x_label: str, x_values, rows +) -> None: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + full_average = [row["full_average_iterations"] for row in rows] + updates_average = [row["update_average_iterations"] for row in rows] + + fig, ax = plt.subplots(figsize=(8, 5)) + ax.plot(x_values, full_average, marker="o", label="Full Recomputation") + ax.plot(x_values, updates_average, marker="o", label="Warm Started Balls") + ax.set_title(title) + ax.set_xlabel(x_label) + ax.set_ylabel("Avg. # of Iterations") + ax.grid(True, which="both", alpha=0.3) + ax.legend() + fig.tight_layout() + fig.savefig(path, dpi=200) + plt.close(fig) + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("--output-dir", default="benchmark_results") + parser.add_argument("--seed", type=int, default=0) + parser.add_argument("--initial-count", type=int, default=100) + parser.add_argument("--distance-factor", type=float, default=1.05) + args = parser.parse_args() + + output_dir = Path(args.output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + rng = np.random.default_rng(args.seed) + + dimension_rows = [ + run_case(rng, dim, 500, args.initial_count, args.distance_factor) + for dim in [5, 10, 50, 100, 500, 1000] + ] + write_csv(output_dir / "vary_dimension.csv", dimension_rows) + print_table("Vary dimension, 500 point insertions", dimension_rows, "dimension") + plot_results( + output_dir / "vary_dimension.png", + "Runtime vs dimension, 500 point insertions", + "# of Dimensions", + [row["dimension"] for row in dimension_rows], + dimension_rows, + ) + plot_iteration_results( + output_dir / "vary_dimension_iterations.png", + "Iterations vs dimension, 500 point insertions", + "# of Dimensions", + [row["dimension"] for row in dimension_rows], + dimension_rows, + ) + + addition_rows = [ + run_case(rng, 500, additions, args.initial_count, args.distance_factor) + for additions in [5, 10, 50, 100, 500, 1000] + ] + write_csv(output_dir / "vary_additions.csv", addition_rows) + print_table("Vary point insertions, 500 dimensions", addition_rows, "additions") + plot_results( + output_dir / "vary_additions.png", + "Runtime vs point insertions, 500 dimensions", + "# of Points Added", + [row["additions"] for row in addition_rows], + addition_rows, + ) + plot_iteration_results( + output_dir / "vary_additions_iterations.png", + "Iterations vs point insertions, 500 dimensions", + "# of Points Added", + [row["additions"] for row in addition_rows], + addition_rows, + ) + + print(f"\nWrote results to {output_dir.resolve()}") + + +if __name__ == "__main__": + main() diff --git a/python/miniball_python.cpp b/python/miniball_python.cpp index 27404a9..1a021c1 100644 --- a/python/miniball_python.cpp +++ b/python/miniball_python.cpp @@ -23,6 +23,7 @@ nb::dict miniball_result(NativeMiniball &mb, size_t dim) { nb::ndarray(mb.center_begin(), {dim}).cast(); result["radius"] = mb.radius(); result["radius_squared"] = mb.squared_radius(); + result["iterations"] = mb.iterations(); return result; } From 222b8bfd59ef63325926e1ece0f91ed95e6d2086 Mon Sep 17 00:00:00 2001 From: Riddho Ridwanul Haque <33618344+RiddhoHaque@users.noreply.github.com> Date: Mon, 11 May 2026 14:10:14 -0400 Subject: [PATCH 3/3] Add SIMD accelerations to QR and Distance Computations --- cpp/main/Seb-inl.h | 98 +++---- cpp/main/Seb.h | 18 +- cpp/main/Seb_point.h | 10 + cpp/main/Seb_simd.h | 452 ++++++++++++++++++++++++++++++++ cpp/main/Subspan-inl.h | 66 ++--- cpp/main/Subspan.h | 4 +- python/benchmark_incremental.py | 148 ++++++++++- python/miniball/__init__.py | 24 +- python/miniball_python.cpp | 28 +- python/test/test_miniball.py | 57 +++- 10 files changed, 770 insertions(+), 135 deletions(-) create mode 100644 cpp/main/Seb_simd.h diff --git a/cpp/main/Seb-inl.h b/cpp/main/Seb-inl.h index 98c74dd..5f6578d 100644 --- a/cpp/main/Seb-inl.h +++ b/cpp/main/Seb-inl.h @@ -54,9 +54,8 @@ namespace SEB_NAMESPACE { unsigned int farthest = 0; // Note: assignment prevents compiler warnings. for (unsigned int j = 1; j < S.size(); ++j) { // compute squared distance from center to S[j]: - Float dist = 0; - for (unsigned int i = 0; i < dim; ++i) - dist += sqr(S[j][i] - center[i]); + Float dist = detail::squared_distance(S[j], center, dim, + use_simd); // enlarge radius if needed: if (dist >= radius_square) { @@ -69,7 +68,8 @@ namespace SEB_NAMESPACE { // initialize support to the farthest point: if (support != NULL) delete support; - support = new Subspan(dim,S,farthest); + support = new Subspan(dim,S,farthest, + use_simd); // statistics: // initialize entry-counters to zero: @@ -110,25 +110,20 @@ namespace SEB_NAMESPACE { Float Smallest_enclosing_ball::find_stop_fraction( const Float* direction, int& stopper) { - using std::inner_product; - Float scale = 0; stopper = -1; const Pt& support_point = S[support->any_member()]; for (unsigned int j = 0; j < S.size(); ++j) if (!support->is_member(j)) { - Float dist = 0; - for (unsigned int i = 0; i < dim; ++i) { - center_to_point[i] = S[j][i] - center[i]; - dist += sqr(center_to_point[i]); - } + Float dist = detail::squared_distance(S[j], center, dim, + use_simd); - for (unsigned int i = 0; i < dim; ++i) - center_to_point[i] = S[j][i] - support_point[i]; + detail::assign_difference(center_to_point, S[j], + support_point, dim, use_simd); const Float denom = - 2 * inner_product(direction, direction+dim, - center_to_point, Float(0)); + 2 * detail::dot(direction, center_to_point, dim, + use_simd); if (denom == 0) continue; @@ -152,8 +147,6 @@ namespace SEB_NAMESPACE { // Further, stopper is set to the index of the most restricting point // and to -1 if no such point was found. { - using std::inner_product; - // We would like to walk the full length of center_to_aff ... Float scale = 1; stopper = -1; @@ -165,12 +158,11 @@ namespace SEB_NAMESPACE { if (!support->is_member(j)) { // compute vector center_to_point from center to the point S[i]: - for (unsigned int i = 0; i < dim; ++i) - center_to_point[i] = S[j][i] - center[i]; + detail::assign_difference(center_to_point, S[j], center, dim, + use_simd); const Float dir_point_prod - = inner_product(center_to_aff,center_to_aff+dim, - center_to_point,Float(0)); + = detail::dot(center_to_aff, center_to_point, dim, use_simd); // we can ignore points beyond support since they stay // enclosed anyway: @@ -184,8 +176,7 @@ namespace SEB_NAMESPACE { // (Better don't try to understand this calculus from the code, // it needs some pencil-and-paper work.) Float bound = radius_square; - bound -= inner_product(center_to_point,center_to_point+dim, - center_to_point,Float(0)); + bound -= detail::squared_norm(center_to_point, dim, use_simd); bound /= 2 * (dist_to_aff_square - dir_point_prod); // watch for numerical instability - if bound=0 then we are @@ -268,14 +259,12 @@ namespace SEB_NAMESPACE { // stopping point exists // walk as far as we can - for (unsigned int i = 0; i < dim; ++i) - center[i] += scale * center_to_aff[i]; + detail::axpy_inplace(center, scale, center_to_aff, dim, use_simd); // update the radius const Pt& stop_point = S[support->any_member()]; - radius_square = 0; - for (unsigned int i = 0; i < dim; ++i) - radius_square += sqr(stop_point[i] - center[i]); + radius_square = detail::squared_distance(stop_point, center, + dim, use_simd); radius_ = sqrt(radius_square); SEB_LOG ("debug"," current radius = " << std::setiosflags(std::ios::scientific) @@ -290,14 +279,12 @@ namespace SEB_NAMESPACE { else { // we can run unhindered into the affine hull SEB_LOG ("debug"," moving into affine hull" << std::endl); - for (unsigned int i=0; iany_member()]; - radius_square = 0; - for (unsigned int i = 0; i < dim; ++i) - radius_square += sqr(stop_point[i] - center[i]); + radius_square = detail::squared_distance(stop_point, center, + dim, use_simd); radius_ = sqrt(radius_square); SEB_LOG ("debug"," current radius = " << std::setiosflags(std::ios::scientific) @@ -336,21 +323,23 @@ namespace SEB_NAMESPACE { SEB_ASSERT(S.size() > 0); SEB_ASSERT(new_point_index < S.size()); const Pt& new_point = S[new_point_index]; - Float dist_square = 0; - for (unsigned int i = 0; i < dim; ++i) - dist_square += sqr(new_point[i] - previous_center[i]); + Float dist_square = detail::squared_distance(new_point, + previous_center, + dim, use_simd); const Float previous_radius = sqrt(previous_squared_radius); const Float dist = sqrt(dist_square); const Float new_radius = (dist + previous_radius) / 2; const Float shift = (dist - previous_radius) / (2 * dist); - for (unsigned int i = 0; i < dim; ++i) - center[i] = previous_center[i] + shift * (new_point[i] - previous_center[i]); + detail::blend_from(center, previous_center, new_point, shift, + dim, use_simd); radius_square = sqr(new_radius); radius_ = new_radius; if (support != NULL) support->reset(new_point_index); else - support = new Subspan(dim, S, new_point_index); + support = new Subspan(dim, S, + new_point_index, + use_simd); SEB_STATS(entry_count = std::vector(S.size(),0)); } else { init_ball(); @@ -376,9 +365,8 @@ namespace SEB_NAMESPACE { if (support != NULL) support->resize_membership(); - Float dist = 0; - for (unsigned int i = 0; i < dim; ++i) - dist += sqr(S[new_point_index][i] - center[i]); + Float dist = detail::squared_distance(S[new_point_index], center, + dim, use_simd); if (dist <= radius_square) return; @@ -392,11 +380,10 @@ namespace SEB_NAMESPACE { support->shortest_vector_to_span(center, center_to_aff); support->shortest_vector_to_span(S[new_point_index], center_to_point); - dist_to_aff_square = 0; - for (unsigned int i = 0; i < dim; ++i) { - center_to_aff[i] -= center_to_point[i]; - dist_to_aff_square += sqr(center_to_aff[i]); - } + dist_to_aff_square = + detail::subtract_inplace_and_squared_norm(center_to_aff, + center_to_point, dim, + use_simd); dist_to_aff = sqrt(dist_to_aff_square); if (dist_to_aff <= Eps * radius_) { update(); @@ -410,13 +397,11 @@ namespace SEB_NAMESPACE { return; } - for (unsigned int i = 0; i < dim; ++i) - center[i] += scale * center_to_aff[i]; + detail::axpy_inplace(center, scale, center_to_aff, dim, use_simd); const Pt& stop_point = S[support->any_member()]; - radius_square = 0; - for (unsigned int i = 0; i < dim; ++i) - radius_square += sqr(stop_point[i] - center[i]); + radius_square = detail::squared_distance(stop_point, center, + dim, use_simd); radius_ = sqrt(radius_square); support->add_point(stopper); @@ -428,7 +413,6 @@ namespace SEB_NAMESPACE { template void Smallest_enclosing_ball::verify() { - using std::inner_product; using std::abs; using std::cout; using std::endl; @@ -449,10 +433,10 @@ namespace SEB_NAMESPACE { for (unsigned int k = 0; k < S.size(); ++k) { // compare center-to-point distance with radius - for (unsigned int i = 0; i < dim; ++i) - center_to_point[i] = S[k][i] - center[i]; - ball_error = sqrt(inner_product(center_to_point,center_to_point+dim, - center_to_point,Float(0))) + detail::assign_difference(center_to_point, S[k], center, dim, + use_simd); + ball_error = sqrt(detail::squared_norm(center_to_point, dim, + use_simd)) - radius_; // check for sphere violations diff --git a/cpp/main/Seb.h b/cpp/main/Seb.h index db8241e..4f6607f 100644 --- a/cpp/main/Seb.h +++ b/cpp/main/Seb.h @@ -9,6 +9,7 @@ #include #include "Seb_configure.h" #include "Seb_point.h" +#include "Seb_simd.h" #include "Subspan.h" namespace SEB_NAMESPACE { @@ -31,11 +32,13 @@ namespace SEB_NAMESPACE { public: // construction and destruction: - Smallest_enclosing_ball(unsigned int d, const PointAccessor &P) + Smallest_enclosing_ball(unsigned int d, const PointAccessor &P, + bool use_simd_if_available = true) // Constructs an instance representing the miniball of points from // set S. The dimension of the ambient space is fixed to d for // lifetime of the instance. : dim(d), S(P), up_to_date(true), support(NULL), + use_simd(use_simd_if_available && detail::simd_available_for()), last_iteration_count(0) { allocate_resources(); @@ -46,11 +49,13 @@ namespace SEB_NAMESPACE { Smallest_enclosing_ball(unsigned int d, const PointAccessor &P, const Float* previous_center, Float previous_squared_radius, - unsigned int new_point_index) + unsigned int new_point_index, + bool use_simd_if_available = true) // Constructs an instance representing the miniball of points from // set S, using the miniball of S without new_point_index as a warm // start. The new point is assumed to lie outside that previous ball. : dim(d), S(P), up_to_date(true), support(NULL), + use_simd(use_simd_if_available && detail::simd_available_for()), last_iteration_count(0) { allocate_resources(); @@ -148,11 +153,15 @@ namespace SEB_NAMESPACE { SEB_ASSERT(!is_empty()); Float dist = 0; - for (unsigned int i = 0; i < dim; ++i) - dist += sqr(point[i] - center[i]); + dist = detail::squared_distance(point, center, dim, use_simd); return dist <= radius_square; } + static bool simd_available() + { + return detail::simd_available(); + } + unsigned int iterations() // Returns the number of iterations used by the most recent update. { @@ -204,6 +213,7 @@ namespace SEB_NAMESPACE { Subspan *support; // the points that lie on the current // boundary and "support" the ball; // the essential structure for update() + bool use_simd; private: // member fields for temporary use: Float *center_to_aff; diff --git a/cpp/main/Seb_point.h b/cpp/main/Seb_point.h index 3896c10..22091e0 100644 --- a/cpp/main/Seb_point.h +++ b/cpp/main/Seb_point.h @@ -63,6 +63,16 @@ namespace SEB_NAMESPACE { return c.end(); } + const Float* data() const + { + return c.empty() ? 0 : &c[0]; + } + + Float* data() + { + return c.empty() ? 0 : &c[0]; + } + private: // member fields: std::vector c; // Cartesian center coordinates }; diff --git a/cpp/main/Seb_simd.h b/cpp/main/Seb_simd.h new file mode 100644 index 0000000..5139abc --- /dev/null +++ b/cpp/main/Seb_simd.h @@ -0,0 +1,452 @@ +// Synopsis: Small vector kernels with optional SIMD acceleration. + +#ifndef SEB_SIMD_H +#define SEB_SIMD_H + +#include +#include "Seb_point.h" + +#if defined(__SSE2__) || defined(_M_X64) || \ + (defined(_M_IX86_FP) && _M_IX86_FP >= 2) +#define SEB_HAS_SSE2 1 +#include +#if defined(_MSC_VER) +#include +#endif +#else +#define SEB_HAS_SSE2 0 +#endif + +namespace SEB_NAMESPACE { +namespace detail { + + template + inline Float local_sqr(const Float x) + { + return x * x; + } + + inline bool simd_available() + { +#if SEB_HAS_SSE2 + // SSE2 is part of the x86-64 ABI. +#if defined(__x86_64__) || defined(_M_X64) + return true; +#elif defined(_MSC_VER) + int info[4]; + __cpuid(info, 1); + return (info[3] & (1 << 26)) != 0; +#elif defined(__GNUC__) || defined(__clang__) + return __builtin_cpu_supports("sse2"); +#else + return true; +#endif +#else + return false; +#endif + } + + template + inline bool simd_available_for() + { + return false; + } + + template<> + inline bool simd_available_for() + { + return simd_available(); + } + + template + inline const Float* contiguous_data(const T&) + { + return 0; + } + + template + inline const Float* contiguous_data(const Float* p) + { + return p; + } + + template + inline const Float* contiguous_data(Float* p) + { + return p; + } + + template + inline const Float* contiguous_data(const Point& p) + { + return p.data(); + } + + template + struct VectorOps { + static Float dot(const Float* a, const Float* b, std::size_t n, bool) + { + Float sum = 0; + for (std::size_t i = 0; i < n; ++i) + sum += a[i] * b[i]; + return sum; + } + + static Float squared_norm(const Float* a, std::size_t n, bool use_simd) + { + return dot(a, a, n, use_simd); + } + + static Float squared_distance(const Float* a, const Float* b, + std::size_t n, bool) + { + Float sum = 0; + for (std::size_t i = 0; i < n; ++i) + sum += local_sqr(a[i] - b[i]); + return sum; + } + + static void assign_difference(Float* out, const Float* a, const Float* b, + std::size_t n, bool) + { + for (std::size_t i = 0; i < n; ++i) + out[i] = a[i] - b[i]; + } + + static void axpy_inplace(Float* y, Float alpha, const Float* x, + std::size_t n, bool) + { + for (std::size_t i = 0; i < n; ++i) + y[i] += alpha * x[i]; + } + + static Float subtract_inplace_and_squared_norm(Float* y, const Float* x, + std::size_t n, bool) + { + Float sum = 0; + for (std::size_t i = 0; i < n; ++i) { + y[i] -= x[i]; + sum += local_sqr(y[i]); + } + return sum; + } + + static void blend(Float* out, const Float* a, const Float* b, Float alpha, + std::size_t n, bool) + { + for (std::size_t i = 0; i < n; ++i) + out[i] = a[i] + alpha * (b[i] - a[i]); + } + + static void rotate_pair(Float* a, Float* b, Float c, Float s, + std::size_t n, bool) + { + for (std::size_t i = 0; i < n; ++i) { + const Float ai = a[i]; + const Float bi = b[i]; + a[i] = c * ai + s * bi; + b[i] = c * bi - s * ai; + } + } + }; + + template<> + struct VectorOps { + static double dot(const double* a, const double* b, std::size_t n, + bool use_simd) + { +#if SEB_HAS_SSE2 + if (use_simd && n >= 2) { + std::size_t i = 0; + __m128d acc = _mm_setzero_pd(); + for (; i + 1 < n; i += 2) + acc = _mm_add_pd(acc, _mm_mul_pd(_mm_loadu_pd(a + i), + _mm_loadu_pd(b + i))); + double tmp[2]; + _mm_storeu_pd(tmp, acc); + double sum = tmp[0] + tmp[1]; + for (; i < n; ++i) + sum += a[i] * b[i]; + return sum; + } +#endif + double sum = 0; + for (std::size_t i = 0; i < n; ++i) + sum += a[i] * b[i]; + return sum; + } + + static double squared_norm(const double* a, std::size_t n, bool use_simd) + { + return dot(a, a, n, use_simd); + } + + static double squared_distance(const double* a, const double* b, + std::size_t n, bool use_simd) + { +#if SEB_HAS_SSE2 + if (use_simd && n >= 2) { + std::size_t i = 0; + __m128d acc = _mm_setzero_pd(); + for (; i + 1 < n; i += 2) { + const __m128d diff = _mm_sub_pd(_mm_loadu_pd(a + i), + _mm_loadu_pd(b + i)); + acc = _mm_add_pd(acc, _mm_mul_pd(diff, diff)); + } + double tmp[2]; + _mm_storeu_pd(tmp, acc); + double sum = tmp[0] + tmp[1]; + for (; i < n; ++i) + sum += local_sqr(a[i] - b[i]); + return sum; + } +#endif + double sum = 0; + for (std::size_t i = 0; i < n; ++i) + sum += local_sqr(a[i] - b[i]); + return sum; + } + + static void assign_difference(double* out, const double* a, + const double* b, std::size_t n, + bool use_simd) + { +#if SEB_HAS_SSE2 + if (use_simd && n >= 2) { + std::size_t i = 0; + for (; i + 1 < n; i += 2) + _mm_storeu_pd(out + i, _mm_sub_pd(_mm_loadu_pd(a + i), + _mm_loadu_pd(b + i))); + for (; i < n; ++i) + out[i] = a[i] - b[i]; + return; + } +#endif + for (std::size_t i = 0; i < n; ++i) + out[i] = a[i] - b[i]; + } + + static void axpy_inplace(double* y, double alpha, const double* x, + std::size_t n, bool use_simd) + { +#if SEB_HAS_SSE2 + if (use_simd && n >= 2) { + std::size_t i = 0; + const __m128d alpha_v = _mm_set1_pd(alpha); + for (; i + 1 < n; i += 2) + _mm_storeu_pd(y + i, + _mm_add_pd(_mm_loadu_pd(y + i), + _mm_mul_pd(alpha_v, _mm_loadu_pd(x + i)))); + for (; i < n; ++i) + y[i] += alpha * x[i]; + return; + } +#endif + for (std::size_t i = 0; i < n; ++i) + y[i] += alpha * x[i]; + } + + static double subtract_inplace_and_squared_norm(double* y, const double* x, + std::size_t n, + bool use_simd) + { +#if SEB_HAS_SSE2 + if (use_simd && n >= 2) { + std::size_t i = 0; + __m128d acc = _mm_setzero_pd(); + for (; i + 1 < n; i += 2) { + const __m128d diff = _mm_sub_pd(_mm_loadu_pd(y + i), + _mm_loadu_pd(x + i)); + _mm_storeu_pd(y + i, diff); + acc = _mm_add_pd(acc, _mm_mul_pd(diff, diff)); + } + double tmp[2]; + _mm_storeu_pd(tmp, acc); + double sum = tmp[0] + tmp[1]; + for (; i < n; ++i) { + y[i] -= x[i]; + sum += local_sqr(y[i]); + } + return sum; + } +#endif + double sum = 0; + for (std::size_t i = 0; i < n; ++i) { + y[i] -= x[i]; + sum += local_sqr(y[i]); + } + return sum; + } + + static void blend(double* out, const double* a, const double* b, + double alpha, std::size_t n, bool use_simd) + { +#if SEB_HAS_SSE2 + if (use_simd && n >= 2) { + std::size_t i = 0; + const __m128d alpha_v = _mm_set1_pd(alpha); + for (; i + 1 < n; i += 2) { + const __m128d av = _mm_loadu_pd(a + i); + const __m128d bv = _mm_loadu_pd(b + i); + _mm_storeu_pd(out + i, + _mm_add_pd(av, _mm_mul_pd(alpha_v, + _mm_sub_pd(bv, av)))); + } + for (; i < n; ++i) + out[i] = a[i] + alpha * (b[i] - a[i]); + return; + } +#endif + for (std::size_t i = 0; i < n; ++i) + out[i] = a[i] + alpha * (b[i] - a[i]); + } + + static void rotate_pair(double* a, double* b, double c, double s, + std::size_t n, bool use_simd) + { +#if SEB_HAS_SSE2 + if (use_simd && n >= 2) { + std::size_t i = 0; + const __m128d c_v = _mm_set1_pd(c); + const __m128d s_v = _mm_set1_pd(s); + for (; i + 1 < n; i += 2) { + const __m128d av = _mm_loadu_pd(a + i); + const __m128d bv = _mm_loadu_pd(b + i); + _mm_storeu_pd(a + i, + _mm_add_pd(_mm_mul_pd(c_v, av), + _mm_mul_pd(s_v, bv))); + _mm_storeu_pd(b + i, + _mm_sub_pd(_mm_mul_pd(c_v, bv), + _mm_mul_pd(s_v, av))); + } + for (; i < n; ++i) { + const double ai = a[i]; + const double bi = b[i]; + a[i] = c * ai + s * bi; + b[i] = c * bi - s * ai; + } + return; + } +#endif + for (std::size_t i = 0; i < n; ++i) { + const double ai = a[i]; + const double bi = b[i]; + a[i] = c * ai + s * bi; + b[i] = c * bi - s * ai; + } + } + }; + + template + inline Float dot(const A& a, const B& b, std::size_t n, bool use_simd) + { + const Float* ap = contiguous_data(a); + const Float* bp = contiguous_data(b); + if (ap != 0 && bp != 0) + return VectorOps::dot(ap, bp, n, + use_simd && simd_available_for()); + + Float sum = 0; + for (std::size_t i = 0; i < n; ++i) + sum += a[i] * b[i]; + return sum; + } + + template + inline Float squared_norm(const A& a, std::size_t n, bool use_simd) + { + const Float* ap = contiguous_data(a); + if (ap != 0) + return VectorOps::squared_norm( + ap, n, use_simd && simd_available_for()); + + Float sum = 0; + for (std::size_t i = 0; i < n; ++i) + sum += local_sqr(a[i]); + return sum; + } + + template + inline Float squared_distance(const A& a, const B& b, std::size_t n, + bool use_simd) + { + const Float* ap = contiguous_data(a); + const Float* bp = contiguous_data(b); + if (ap != 0 && bp != 0) + return VectorOps::squared_distance(ap, bp, n, + use_simd && + simd_available_for()); + + Float sum = 0; + for (std::size_t i = 0; i < n; ++i) + sum += local_sqr(a[i] - b[i]); + return sum; + } + + template + inline void assign_difference(Float* out, const A& a, const B& b, + std::size_t n, bool use_simd) + { + const Float* ap = contiguous_data(a); + const Float* bp = contiguous_data(b); + if (ap != 0 && bp != 0) { + VectorOps::assign_difference(out, ap, bp, n, + use_simd && + simd_available_for()); + return; + } + + for (std::size_t i = 0; i < n; ++i) + out[i] = a[i] - b[i]; + } + + template + inline void axpy_inplace(Float* y, Float alpha, const Float* x, + std::size_t n, bool use_simd) + { + VectorOps::axpy_inplace(y, alpha, x, n, + use_simd && simd_available_for()); + } + + template + inline Float subtract_inplace_and_squared_norm(Float* y, const Float* x, + std::size_t n, bool use_simd) + { + return VectorOps::subtract_inplace_and_squared_norm( + y, x, n, use_simd && simd_available_for()); + } + + template + inline void blend(Float* out, const Float* a, const Float* b, Float alpha, + std::size_t n, bool use_simd) + { + VectorOps::blend(out, a, b, alpha, n, + use_simd && simd_available_for()); + } + + template + inline void blend_from(Float* out, const Float* a, const B& b, Float alpha, + std::size_t n, bool use_simd) + { + const Float* bp = contiguous_data(b); + if (bp != 0) { + blend(out, a, bp, alpha, n, use_simd); + return; + } + + for (std::size_t i = 0; i < n; ++i) + out[i] = a[i] + alpha * (b[i] - a[i]); + } + + template + inline void rotate_pair(Float* a, Float* b, Float c, Float s, + std::size_t n, bool use_simd) + { + VectorOps::rotate_pair(a, b, c, s, n, + use_simd && simd_available_for()); + } + +} // namespace detail +} // namespace SEB_NAMESPACE + +#endif // SEB_SIMD_H diff --git a/cpp/main/Subspan-inl.h b/cpp/main/Subspan-inl.h index d76fd89..87f96b0 100644 --- a/cpp/main/Subspan-inl.h +++ b/cpp/main/Subspan-inl.h @@ -10,6 +10,7 @@ #include #include #include "Seb_configure.h" +#include "Seb_simd.h" #include "Subspan.h" // Note: header included for better syntax highlighting in some IDEs. @@ -51,8 +52,12 @@ namespace SEB_NAMESPACE { } template - Subspan::Subspan(unsigned int dim, const PointAccessor& S, int index) - : S(S), membership(S.size()), dim(dim), members(dim+1) + Subspan::Subspan(unsigned int dim, + const PointAccessor& S, + int index, + bool use_simd_if_available) + : S(S), membership(S.size()), dim(dim), members(dim+1), + use_simd(use_simd_if_available && detail::simd_available_for()) { // allocate storage for Q, R, u, and w: Q = new Float *[dim]; @@ -174,20 +179,16 @@ namespace SEB_NAMESPACE { shortest_vector_to_span(RandomAccessIterator1 p, RandomAccessIterator2 w) { - using std::inner_product; - // compute vector from p to origin, i.e., w = origin - p: - for (unsigned int i=0; i(w, SEB_AFFINE_ORIGIN, p, dim, use_simd); // remove projections of w onto the affine hull: for (unsigned int j = 0; j < r; ++j) { - const Float scale = inner_product(w,w+dim,Q[j],Float(0)); - for (unsigned int i = 0; i < dim; ++i) - w[i] -= scale * Q[j][i]; + const Float scale = detail::dot(w, Q[j], dim, use_simd); + detail::axpy_inplace(w, -scale, Q[j], dim, use_simd); } - return inner_product(w,w+dim,w,Float(0)); + return detail::squared_norm(w, dim, use_simd); } template @@ -230,15 +231,11 @@ namespace SEB_NAMESPACE { RandomAccessIterator2 lambdas) { // compute relative position of p, i.e., u = p - origin: - for (unsigned int i=0; i(u, p, SEB_AFFINE_ORIGIN, dim, use_simd); // calculate Q^T u into w: - for (unsigned int i = 0; i < dim; ++i) { - w[i] = 0; - for (unsigned int k = 0; k < dim; ++k) - w[i] += Q[i][k] * u[k]; - } + for (unsigned int i = 0; i < dim; ++i) + w[i] = detail::dot(Q[i], u, dim, use_simd); // We compute the coefficients by backsubstitution. Notice that // @@ -269,11 +266,8 @@ namespace SEB_NAMESPACE { SEB_ASSERT(r < dim); // compute new column R[r] = Q^T * u - for (unsigned int i = 0; i < dim; ++i) { - R[r][i] = 0; - for (unsigned int k = 0; k < dim; ++k) - R[r][i] += Q[i][k] * u[k]; - } + for (unsigned int i = 0; i < dim; ++i) + R[r][i] = detail::dot(Q[i], u, dim, use_simd); // zero all entries R[r][dim-1] down to R[r][r+1] for (unsigned int j = dim-1; j > r; --j) { @@ -288,12 +282,7 @@ namespace SEB_NAMESPACE { R[r][j-1] = c * R[r][j-1] + s * R[r][j]; // rotate two Q-columns - for (unsigned int i = 0; i < dim; ++i) { - const Float a = Q[j-1][i]; - const Float b = Q[j][i]; - Q[j-1][i] = c * a + s * b; - Q[j][i] = c * b - s * a; - } + detail::rotate_pair(Q[j-1], Q[j], c, s, dim, use_simd); } } @@ -323,12 +312,7 @@ namespace SEB_NAMESPACE { } // rotate Q-columns - for (unsigned int i = 0; i < dim; ++i) { - const Float a = Q[pos][i]; - const Float b = Q[pos+1][i]; - Q[pos][i] = c * a + s * b; - Q[pos+1][i] = c * b - s * a; - } + detail::rotate_pair(Q[pos], Q[pos+1], c, s, dim, use_simd); } } @@ -338,11 +322,8 @@ namespace SEB_NAMESPACE { // A + u * [1,...,1] = Q' R'. { // compute w = Q^T * u - for (unsigned int i = 0; i < dim; ++i) { - w[i] = 0; - for (unsigned int k = 0; k < dim; ++k) - w[i] += Q[i][k] * u[k]; - } + for (unsigned int i = 0; i < dim; ++i) + w[i] = detail::dot(Q[i], u, dim, use_simd); // rotate w down to a multiple of the first unit vector; // the operations have to be recorded in R and Q @@ -370,12 +351,7 @@ namespace SEB_NAMESPACE { } // rotate two Q-columns - for (unsigned int i = 0; i < dim; ++i) { - const Float a = Q[k-1][i]; - const Float b = Q[k][i]; - Q[k-1][i] = c * a + s * b; - Q[k][i] = c * b - s * a; - } + detail::rotate_pair(Q[k-1], Q[k], c, s, dim, use_simd); } // add w * (1,...,1)^T to new R diff --git a/cpp/main/Subspan.h b/cpp/main/Subspan.h index 5690161..dec1f9f 100644 --- a/cpp/main/Subspan.h +++ b/cpp/main/Subspan.h @@ -76,7 +76,8 @@ namespace SEB_NAMESPACE { { public: // construction and deletion: - Subspan(unsigned int dim, const PointAccessor& S, int i); + Subspan(unsigned int dim, const PointAccessor& S, int i, + bool use_simd_if_available = true); // Constructs an instance representing the affine hull aff(M) of M={p}, // where p is the point S[i] from S. // @@ -172,6 +173,7 @@ namespace SEB_NAMESPACE { // in row i and column j Float *u,*w; // needed for rank-1 update unsigned int r; // the rank of R (i.e. #points - 1) + bool use_simd; // use vector kernels when available }; } // namespace SEB_NAMESPACE diff --git a/python/benchmark_incremental.py b/python/benchmark_incremental.py index db40abd..5681822 100644 --- a/python/benchmark_incremental.py +++ b/python/benchmark_incremental.py @@ -8,7 +8,7 @@ import numpy as np -from miniball import Miniball, miniball +from miniball import Miniball, miniball, simd_available def unit_direction(rng: np.random.Generator, dim: int) -> np.ndarray: @@ -27,7 +27,7 @@ def generate_outside_sequence( distance_factor: float, ) -> np.ndarray: points = np.asarray(initial_points, dtype=np.float64) - stream = Miniball(points) + stream = Miniball(points, use_simd_if_available=False) current = stream.result() additions = [] @@ -45,13 +45,15 @@ def generate_outside_sequence( def time_full_recomputation( - all_points: np.ndarray, initial_count: int + all_points: np.ndarray, initial_count: int, use_simd_if_available: bool ) -> tuple[float, list[np.ndarray], dict[str, float]]: centers = [] iterations = [] start = time.perf_counter() for size in range(initial_count + 1, len(all_points) + 1): - current = miniball(all_points[:size]) + current = miniball( + all_points[:size], use_simd_if_available=use_simd_if_available + ) centers.append(np.asarray(current["center"])) iterations.append(int(current["iterations"])) elapsed = time.perf_counter() - start @@ -59,9 +61,14 @@ def time_full_recomputation( def time_incremental( - all_points: np.ndarray, initial_count: int, centers: list[np.ndarray] + all_points: np.ndarray, + initial_count: int, + centers: list[np.ndarray], + use_simd_if_available: bool, ) -> tuple[float, dict[str, float]]: - stream = Miniball(all_points[:initial_count]) + stream = Miniball( + all_points[:initial_count], use_simd_if_available=use_simd_if_available + ) iterations = [] start = time.perf_counter() @@ -96,6 +103,7 @@ def run_case( additions_count: int, initial_count: int, distance_factor: float, + run_simd: bool, ) -> dict[str, float]: initial_points = rng.normal(size=(initial_count, dim)) additions = generate_outside_sequence( @@ -104,7 +112,7 @@ def run_case( all_points = np.vstack([initial_points, additions]) full_seconds, centers, full_iteration_stats = time_full_recomputation( - all_points, initial_count + all_points, initial_count, use_simd_if_available=False ) print('Finished full recomputation for dimension {}, point insertions {}'.format(dim, additions_count)) print('Full recomputation took {:.4f} seconds'.format(full_seconds)) @@ -116,7 +124,7 @@ def run_case( ) ) update_seconds, iteration_stats = time_incremental( - all_points, initial_count, centers + all_points, initial_count, centers, use_simd_if_available=False ) print('Finished updates for dimension {}, point insertions {}'.format(dim, additions_count)) print('Updates took {:.4f} seconds'.format(update_seconds)) @@ -127,7 +135,7 @@ def run_case( iteration_stats["update_average_iterations"], ) ) - return { + row = { "dimension": dim, "additions": additions_count, "full_seconds": full_seconds, @@ -139,6 +147,78 @@ def run_case( **iteration_stats, } + if run_simd: + simd_full_seconds, simd_centers, simd_full_iteration_stats = ( + time_full_recomputation( + all_points, initial_count, use_simd_if_available=True + ) + ) + print( + "Finished SIMD full recomputation for dimension {}, point insertions {}".format( + dim, additions_count + ) + ) + print("SIMD full recomputation took {:.4f} seconds".format(simd_full_seconds)) + print( + "SIMD full recomputation iterations: max {}, median {:.1f}, average {:.2f}".format( + simd_full_iteration_stats["full_max_iterations"], + simd_full_iteration_stats["full_median_iterations"], + simd_full_iteration_stats["full_average_iterations"], + ) + ) + + simd_update_seconds, simd_iteration_stats = time_incremental( + all_points, + initial_count, + simd_centers, + use_simd_if_available=True, + ) + print( + "Finished SIMD updates for dimension {}, point insertions {}".format( + dim, additions_count + ) + ) + print("SIMD updates took {:.4f} seconds".format(simd_update_seconds)) + print( + "SIMD iterations per update: max {}, median {:.1f}, average {:.2f}".format( + simd_iteration_stats["update_max_iterations"], + simd_iteration_stats["update_median_iterations"], + simd_iteration_stats["update_average_iterations"], + ) + ) + + row.update( + { + "simd_full_seconds": simd_full_seconds, + "simd_update_seconds": simd_update_seconds, + "simd_speedup": ( + simd_full_seconds / simd_update_seconds + if simd_update_seconds > 0 + else float("inf") + ), + "simd_full_max_iterations": simd_full_iteration_stats[ + "full_max_iterations" + ], + "simd_full_median_iterations": simd_full_iteration_stats[ + "full_median_iterations" + ], + "simd_full_average_iterations": simd_full_iteration_stats[ + "full_average_iterations" + ], + "simd_update_max_iterations": simd_iteration_stats[ + "update_max_iterations" + ], + "simd_update_median_iterations": simd_iteration_stats[ + "update_median_iterations" + ], + "simd_update_average_iterations": simd_iteration_stats[ + "update_average_iterations" + ], + } + ) + + return row + def write_csv(path: Path, rows: list[dict[str, float]]) -> None: with path.open("w", newline="") as file: @@ -150,13 +230,20 @@ def write_csv(path: Path, rows: list[dict[str, float]]) -> None: def print_table(title: str, rows: list[dict[str, float]], x_key: str) -> None: print(f"\n{title}") print("-" * len(title)) - print( + has_simd = bool(rows and "simd_full_seconds" in rows[0]) + header = ( f"{x_key:>10} {'full (s)':>12} {'updates (s)':>17} " f"{'speedup':>10} {'full max':>8} {'full med':>8} " f"{'full avg':>8} {'upd max':>8} {'upd med':>8} {'upd avg':>8}" ) + if has_simd: + header += ( + f" {'simd full (s)':>14} {'simd updates (s)':>17} " + f"{'simd speedup':>12}" + ) + print(header) for row in rows: - print( + line = ( f"{int(row[x_key]):>10} " f"{row['full_seconds']:>12.4f} " f"{row['update_seconds']:>17.4f} " @@ -168,6 +255,13 @@ def print_table(title: str, rows: list[dict[str, float]], x_key: str) -> None: f"{row['update_median_iterations']:>8.1f} " f"{row['update_average_iterations']:>8.2f}" ) + if has_simd: + line += ( + f" {row['simd_full_seconds']:>14.4f} " + f"{row['simd_update_seconds']:>17.4f} " + f"{row['simd_speedup']:>12.2f}" + ) + print(line) def plot_results(path: Path, title: str, x_label: str, x_values, rows) -> None: @@ -178,10 +272,33 @@ def plot_results(path: Path, title: str, x_label: str, x_values, rows) -> None: full = [row["full_seconds"] for row in rows] updates = [row["update_seconds"] for row in rows] + simd_full = ( + [row["simd_full_seconds"] for row in rows] + if rows and "simd_full_seconds" in rows[0] + else None + ) + simd_updates = ( + [row["simd_update_seconds"] for row in rows] + if rows and "simd_update_seconds" in rows[0] + else None + ) fig, ax = plt.subplots(figsize=(8, 5)) ax.plot(x_values, full, marker="o", label="Full Recomputation") ax.plot(x_values, updates, marker="o", label="Warm Started Balls") + if simd_full is not None and simd_updates is not None: + ax.plot( + x_values, + simd_full, + marker="o", + label="Full Recomputation w/ SIMD acceleration", + ) + ax.plot( + x_values, + simd_updates, + marker="o", + label="Warm Started Balls w/ SIMD acceleration", + ) ax.set_title(title) ax.set_xlabel(x_label) ax.set_ylabel("Runtime (seconds)") @@ -227,9 +344,14 @@ def main() -> None: output_dir = Path(args.output_dir) output_dir.mkdir(parents=True, exist_ok=True) rng = np.random.default_rng(args.seed) + run_simd = simd_available() + if run_simd: + print("SIMD acceleration is available; benchmark will include SIMD runs.") + else: + print("SIMD acceleration is not available; benchmark will run scalar only.") dimension_rows = [ - run_case(rng, dim, 500, args.initial_count, args.distance_factor) + run_case(rng, dim, 500, args.initial_count, args.distance_factor, run_simd) for dim in [5, 10, 50, 100, 500, 1000] ] write_csv(output_dir / "vary_dimension.csv", dimension_rows) @@ -250,7 +372,7 @@ def main() -> None: ) addition_rows = [ - run_case(rng, 500, additions, args.initial_count, args.distance_factor) + run_case(rng, 500, additions, args.initial_count, args.distance_factor, run_simd) for additions in [5, 10, 50, 100, 500, 1000] ] write_csv(output_dir / "vary_additions.csv", addition_rows) diff --git a/python/miniball/__init__.py b/python/miniball/__init__.py index 2af671a..8e07e2f 100644 --- a/python/miniball/__init__.py +++ b/python/miniball/__init__.py @@ -2,14 +2,21 @@ _Miniball, _compute_miniball, _compute_miniball_incremental, + _simd_available, ) import numpy as np -__all__ = ["miniball", "incremental_miniball", "Miniball"] +__all__ = ["miniball", "incremental_miniball", "Miniball", "simd_available"] -def miniball(points: np.typing.ArrayLike): +def simd_available() -> bool: + """Return whether the native extension can use SIMD kernels.""" + + return bool(_simd_available()) + + +def miniball(points: np.typing.ArrayLike, use_simd_if_available: bool = True): """Compute the smallest enclosing ball for a set of points.""" points = np.ascontiguousarray(points, dtype=np.float64) @@ -17,13 +24,14 @@ def miniball(points: np.typing.ArrayLike): if len(points.shape) != 2: msg = f"Input array must be 2-dimensional! Got shape `{points.shape}`." raise TypeError(msg) - return _compute_miniball(points) + return _compute_miniball(points, use_simd_if_available) def incremental_miniball( points: np.typing.ArrayLike, current: dict, point: np.typing.ArrayLike, + use_simd_if_available: bool = True, ): """Compute the miniball after appending one point to an existing point set.""" @@ -47,18 +55,22 @@ def incremental_miniball( return current # Otherwise the native helper starts from a ball with point on the boundary. - return _compute_miniball_incremental(points, center, radius_squared, point) + return _compute_miniball_incremental( + points, center, radius_squared, point, use_simd_if_available + ) class Miniball: """Maintain a miniball as points are appended.""" - def __init__(self, points: np.typing.ArrayLike): + def __init__( + self, points: np.typing.ArrayLike, use_simd_if_available: bool = True + ): points = np.ascontiguousarray(points, dtype=np.float64) if len(points.shape) != 2: msg = f"Input array must be 2-dimensional! Got shape `{points.shape}`." raise TypeError(msg) - self._miniball = _Miniball(points) + self._miniball = _Miniball(points, use_simd_if_available) def add(self, point: np.typing.ArrayLike): point = np.ascontiguousarray(point, dtype=np.float64) diff --git a/python/miniball_python.cpp b/python/miniball_python.cpp index 1a021c1..46f4a0d 100644 --- a/python/miniball_python.cpp +++ b/python/miniball_python.cpp @@ -31,7 +31,9 @@ nb::dict miniball_result(NativeMiniball &mb, size_t dim) { class Miniball { public: Miniball( - nb::ndarray, nb::c_contig> points_arr) { + nb::ndarray, nb::c_contig> points_arr, + bool use_simd_if_available = true) + : use_simd_if_available_(use_simd_if_available) { size_t n_points = points_arr.shape(0); dim_ = points_arr.shape(1); if (n_points == 0) { @@ -44,7 +46,8 @@ class Miniball { for (size_t i = 0; i < n_points; ++i) { points_->emplace_back(dim_, data + i * dim_); } - miniball_.reset(new NativeMiniball(static_cast(dim_), *points_)); + miniball_.reset(new NativeMiniball(static_cast(dim_), *points_, + use_simd_if_available_)); } nb::dict add(nb::ndarray, nb::c_contig> point_arr) { @@ -76,6 +79,7 @@ class Miniball { } size_t dim_; + bool use_simd_if_available_; std::unique_ptr> points_; std::unique_ptr miniball_; }; @@ -92,7 +96,8 @@ class Miniball { * and "radius_squared" (float). */ nb::dict compute_miniball( - nb::ndarray, nb::c_contig> points_arr) { + nb::ndarray, nb::c_contig> points_arr, + bool use_simd_if_available = true) { size_t n_points = points_arr.shape(0); size_t dim = points_arr.shape(1); @@ -107,7 +112,7 @@ nb::dict compute_miniball( } // Compute the smallest enclosing ball. - NativeMiniball mb(dim, points); + NativeMiniball mb(dim, points, use_simd_if_available); return miniball_result(mb, dim); } @@ -123,7 +128,8 @@ nb::dict compute_miniball_incremental( nb::ndarray, nb::c_contig> points_arr, nb::ndarray, nb::c_contig> center_arr, double radius_squared, - nb::ndarray, nb::c_contig> point_arr) { + nb::ndarray, nb::c_contig> point_arr, + bool use_simd_if_available = true) { size_t n_points = points_arr.shape(0); size_t dim = points_arr.shape(1); @@ -144,7 +150,8 @@ nb::dict compute_miniball_incremental( points.emplace_back(dim, new_point); NativeMiniball mb(dim, points, old_center, radius_squared, - static_cast(n_points)); + static_cast(n_points), + use_simd_if_available); return miniball_result(mb, dim); } @@ -153,13 +160,18 @@ nb::dict compute_miniball_incremental( // This replaces all the PyMethodDef, PyModuleDef, and PyInit boilerplate. NB_MODULE(_miniball, m) { m.def("_compute_miniball", &compute_miniball, nb::arg("points"), + nb::arg("use_simd_if_available") = true, "Compute the smallest enclosing ball for a set of points."); m.def("_compute_miniball_incremental", &compute_miniball_incremental, nb::arg("points"), nb::arg("center"), nb::arg("radius_squared"), - nb::arg("point"), + nb::arg("point"), nb::arg("use_simd_if_available") = true, "Compute the smallest enclosing ball after appending one outside point."); + m.def("_simd_available", &NativeMiniball::simd_available, + "Return whether this build can use SIMD kernels on this CPU."); nb::class_(m, "_Miniball") - .def(nb::init, nb::c_contig>>()) + .def(nb::init, nb::c_contig>, + bool>(), + nb::arg("points"), nb::arg("use_simd_if_available") = true) .def("add", &Miniball::add) .def("add_points", &Miniball::add_points) .def("result", &Miniball::result); diff --git a/python/test/test_miniball.py b/python/test/test_miniball.py index b59ae6e..9c46b3e 100644 --- a/python/test/test_miniball.py +++ b/python/test/test_miniball.py @@ -1,5 +1,5 @@ import pytest -from miniball import Miniball, incremental_miniball, miniball +from miniball import Miniball, incremental_miniball, miniball, simd_available import numpy as np @@ -122,3 +122,58 @@ def test_stateful_incremental_dimension_mismatch(): def test_stateful_incremental_rejects_empty_initial_points(): with pytest.raises(Exception, match="at least one point"): Miniball(np.empty((0, 2), dtype=np.double)) + + +def test_simd_availability_flag_is_boolean(): + assert isinstance(simd_available(), bool) + + +def test_simd_and_scalar_find_same_ball_for_random_high_dimensional_points(): + rng = np.random.default_rng(42) + points = rng.normal(size=(300, 64)) + + scalar = miniball(points, use_simd_if_available=False) + simd = miniball(points, use_simd_if_available=True) + + np.testing.assert_allclose(simd["center"], scalar["center"], rtol=1e-10) + assert simd["radius"] == pytest.approx(scalar["radius"], rel=1e-10) + assert simd["radius_squared"] == pytest.approx( + scalar["radius_squared"], rel=1e-10 + ) + + +def test_simd_and_scalar_incremental_find_same_ball(): + rng = np.random.default_rng(43) + points = rng.normal(size=(200, 32)) + point = rng.normal(size=32) * 4 + scalar_current = miniball(points, use_simd_if_available=False) + simd_current = miniball(points, use_simd_if_available=True) + + scalar = incremental_miniball( + points, scalar_current, point, use_simd_if_available=False + ) + simd = incremental_miniball(points, simd_current, point, use_simd_if_available=True) + + np.testing.assert_allclose(simd["center"], scalar["center"], rtol=1e-10) + assert simd["radius"] == pytest.approx(scalar["radius"], rel=1e-10) + assert simd["radius_squared"] == pytest.approx( + scalar["radius_squared"], rel=1e-10 + ) + + +def test_simd_and_scalar_stateful_incremental_find_same_ball(): + rng = np.random.default_rng(44) + initial = rng.normal(size=(100, 24)) + additions = rng.normal(size=(20, 24)) + scalar = Miniball(initial, use_simd_if_available=False) + simd = Miniball(initial, use_simd_if_available=True) + + for point in additions: + scalar_result = scalar.add(point) + simd_result = simd.add(point) + + np.testing.assert_allclose(simd_result["center"], scalar_result["center"], rtol=1e-10) + assert simd_result["radius"] == pytest.approx(scalar_result["radius"], rel=1e-10) + assert simd_result["radius_squared"] == pytest.approx( + scalar_result["radius_squared"], rel=1e-10 + )