From b2e8c07680d2cd85e02aaa4ca27dfc4a8494d5c6 Mon Sep 17 00:00:00 2001 From: ADEGA Date: Tue, 29 Apr 2025 09:28:48 -0500 Subject: [PATCH 1/8] First start low-variance sampling Signed-off-by: ADEGA --- beluga/include/beluga/views.hpp | 1 + .../beluga/views/low_variance_sample.hpp | 120 ++++++++++++++++++ 2 files changed, 121 insertions(+) create mode 100644 beluga/include/beluga/views/low_variance_sample.hpp diff --git a/beluga/include/beluga/views.hpp b/beluga/include/beluga/views.hpp index f530d50a5a..5bf271c925 100644 --- a/beluga/include/beluga/views.hpp +++ b/beluga/include/beluga/views.hpp @@ -16,6 +16,7 @@ #define BELUGA_VIEWS_HPP #include +#include #include #include #include diff --git a/beluga/include/beluga/views/low_variance_sample.hpp b/beluga/include/beluga/views/low_variance_sample.hpp new file mode 100644 index 0000000000..e15732ec10 --- /dev/null +++ b/beluga/include/beluga/views/low_variance_sample.hpp @@ -0,0 +1,120 @@ +// Copyright 2023-2025 Ekumen, Inc. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef BELUGA_VIEWS_LOW_VARIANCE_SAMPLE_HPP +#define BELUGA_VIEWS_LOW_VARIANCE_SAMPLE_HPP + +#include +#include +#include +#include +#include +#include + +#include + +#include + +namespace beluga { + +namespace views { + +namespace detail { + +/// Implementation detail for a low_variance_sample range adaptor object. +struct low_variance_sample_fn { + /// Overload that implements the low_variance_sample algorithm. + /** + * \tparam Range An [input range](https://en.cppreference.com/w/cpp/ranges/input_range) with particle states + * or weights. + * \param range Source range from where to sample elements. + * \param n The number of samples to take (M in the algorithm). + */ + template , int> = 0> + constexpr auto operator()(Range&& range, std::size_t n) const { + static_assert(ranges::input_range); + + using RangeValue = ranges::range_value_t; + std::vector weights; + std::vector particles; // Store particles for direct access + + if constexpr (is_particle_range_v) { + weights.reserve(ranges::distance(range)); + particles.reserve(ranges::distance(range)); + ranges::for_each(range, [&](const RangeValue& particle) { + weights.push_back(beluga::weight(particle)); + particles.push_back(particle); + }); + } else { + weights = ranges::to>(range); + particles = ranges::to>(range); + } + + const std::size_t num_particles = weights.size(); + if (num_particles == 0 || n == 0) { + return ranges::empty(range); + } + + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_real_distribution<> distrib(0.0, 1.0); + const double r = distrib(gen); + const double inv_n = 1.0 / static_cast(n); + double c = weights[0]; + std::size_t i = 0; + std::vector resampled_particles; + resampled_particles.reserve(n); + + for (std::size_t m = 1; m <= n; ++m) { + const double u = r + static_cast(m - 1) * inv_n; + while (u > c) { + i++; + if (i < num_particles) { + c += weights[i]; + } else { + i = num_particles - 1; + break; + } + } + resampled_particles.push_back(particles[i]); + } + + return resampled_particles | ranges::views::common; + } + + /// Overload that returns a view closure to compose with other views. + /** + * \param n The number of samples to take. + */ + constexpr auto operator()(std::size_t n) const { + return ranges::make_view_closure(ranges::bind_back(low_variance_sample_fn{}, n)); + } +}; + +} // namespace detail + +/// [Range adaptor object](https://en.cppreference.com/w/cpp/named_req/RangeAdaptorObject) that +/// will perform low-variance sampling on the input range using the provided algorithm. +/** + * The input range should be a set of particles with associated weights. + * This adaptor will return a new range of resampled particles using the low-variance sampling + * algorithm described. + */ +inline constexpr detail::low_variance_sample_fn low_variance_sample; + +} // namespace views + +} // namespace beluga + +#endif From 99341f62e30a9b944719a79d03bc89834600e493 Mon Sep 17 00:00:00 2001 From: ADEGA Date: Mon, 12 May 2025 14:47:18 -0500 Subject: [PATCH 2/8] [Refactor] Addressing Comments. Making low_variance_sample similar to sample Signed-off-by: ADEGA --- .../beluga/views/low_variance_sample.hpp | 275 +++++++++++++----- 1 file changed, 201 insertions(+), 74 deletions(-) diff --git a/beluga/include/beluga/views/low_variance_sample.hpp b/beluga/include/beluga/views/low_variance_sample.hpp index e15732ec10..99e8105ba7 100644 --- a/beluga/include/beluga/views/low_variance_sample.hpp +++ b/beluga/include/beluga/views/low_variance_sample.hpp @@ -1,4 +1,4 @@ -// Copyright 2023-2025 Ekumen, Inc. +// Copyright 2024 Ekumen, Inc. // // Licensed under the Apache License, Version 2.0 (the "License"); // you may not use this file except in compliance with the License. @@ -15,106 +15,233 @@ #ifndef BELUGA_VIEWS_LOW_VARIANCE_SAMPLE_HPP #define BELUGA_VIEWS_LOW_VARIANCE_SAMPLE_HPP -#include -#include -#include -#include #include -#include -#include +#include +#include +#include #include +#include -namespace beluga { +/** + * \file + * \brief Implementation of a sample (with replacement) range adaptor object. + */ -namespace views { +namespace beluga::views { namespace detail { -/// Implementation detail for a low_variance_sample range adaptor object. -struct low_variance_sample_fn { - /// Overload that implements the low_variance_sample algorithm. +/// Implementation of the sample view. +/** + * \tparam Range A [random access](https://en.cppreference.com/w/cpp/ranges/random_access_range) and + * [sized](https://en.cppreference.com/w/cpp/ranges/sized_range) range. + * \tparam Weights + * \tparam URNG A random number generator that satisfies the + * [UniformRandomBitGenerator](https://en.cppreference.com/w/cpp/named_req/UniformRandomBitGenerator) + * requirements. + */ +template +struct low_variance_sample_view + : public ranges::view_facade, ranges::infinite> { + public: + /// Default constructor. + low_variance_sample_view() = default; + + /// Construct the view from an existing range. /** - * \tparam Range An [input range](https://en.cppreference.com/w/cpp/ranges/input_range) with particle states - * or weights. - * \param range Source range from where to sample elements. - * \param n The number of samples to take (M in the algorithm). + * \param range The range to be adapted. + * \param weights The weights associated with the elements in the range. + * \param engine The random number generator object. */ - template , int> = 0> - constexpr auto operator()(Range&& range, std::size_t n) const { - static_assert(ranges::input_range); - - using RangeValue = ranges::range_value_t; - std::vector weights; - std::vector particles; // Store particles for direct access - - if constexpr (is_particle_range_v) { - weights.reserve(ranges::distance(range)); - particles.reserve(ranges::distance(range)); - ranges::for_each(range, [&](const RangeValue& particle) { - weights.push_back(beluga::weight(particle)); - particles.push_back(particle); - }); - } else { - weights = ranges::to>(range); - particles = ranges::to>(range); - } - - const std::size_t num_particles = weights.size(); - if (num_particles == 0 || n == 0) { - return ranges::empty(range); - } + constexpr low_variance_sample_view(Range range, Weights weights, URNG& engine = ranges::detail::get_random_engine()) + : range_{std::move(range)}, weights_{std::move(weights)}, engine_{std::addressof(engine)} { + assert(ranges::size(range) > 0); + } - std::random_device rd; - std::mt19937 gen(rd()); - std::uniform_real_distribution<> distrib(0.0, 1.0); - const double r = distrib(gen); - const double inv_n = 1.0 / static_cast(n); - double c = weights[0]; - std::size_t i = 0; - std::vector resampled_particles; - resampled_particles.reserve(n); - - for (std::size_t m = 1; m <= n; ++m) { - const double u = r + static_cast(m - 1) * inv_n; - while (u > c) { - i++; - if (i < num_particles) { - c += weights[i]; - } else { - i = num_particles - 1; - break; - } + private: + // `ranges::range_access` needs access to the cursor members. + friend ranges::range_access; + + static_assert(ranges::sized_range); + static_assert(ranges::random_access_range); + static_assert(ranges::input_range); + + /// Cursor class that handles the iteration logic. + struct cursor { + public: + /// Default constructor. + cursor() = default; + + /// Construct a cursor from the parent view elements. + constexpr explicit cursor(low_variance_sample_view* view) + : view_(view), + range_begin_{ranges::begin(view_->range_)}, + weights_begin_{ranges::begin(view_->weights_)}, + M_{static_cast(ranges::size(view_->range_))}, + r_{std::uniform_real_distribution(0.0, 1.0 / static_cast(M_))( + *view_->engine_)}, // Generate r only once + m_{0}, + i_{0}, + c_{*weights_begin_} {} + + /// Access the current iterator. + [[nodiscard]] constexpr decltype(auto) read() const noexcept(noexcept(*range_begin_)) { return *it_; } + + /// Position the current iterator. + constexpr void next() { + ++m_; + const double U = r_ + (static_cast(m_) - 1.0) / static_cast(M_); + while (U > c_) { + ++i_; + c_ += static_cast(*std::next(weights_begin_, i_)); // weights_begin_[i_] } - resampled_particles.push_back(particles[i]); + it_ = std::next(range_begin_, i_); // range_begin_ + i_ } - return resampled_particles | ranges::views::common; + /// Returns true if the cursor is at the end of the range. + constexpr bool equal(const cursor& other) const noexcept { return m_ == other.M_; } + + private: + low_variance_sample_view* view_; + ranges::iterator_t range_begin_; + ranges::iterator_t it_; + ranges::iterator_t weights_begin_; + long long M_; + double r_; + long long m_; + long long i_; + double c_; + }; + + /// Return the cursor for the begin iterator. + [[nodiscard]] constexpr auto begin_cursor() { return cursor{this}; } + + /// Return an unreachable sentinel since this is an infinite range. + [[nodiscard]] constexpr auto end_cursor() const noexcept { return ranges::unreachable_sentinel_t{}; } + + Range range_; + Weights weights_; + URNG* engine_; +}; + +/// Implementation detail for a low_variance_sample algorithm. +struct low_variance_sample_base_fn { + protected: + /// Sample from weighted ranges. + template + constexpr auto low_variance_sample_from_range(Range&& range, Weights&& weights, URNG& engine) const { + static_assert(ranges::sized_range); + static_assert(ranges::random_access_range); + static_assert(ranges::input_range); + return low_variance_sample_view{ranges::views::all(std::forward(range)), std::move(weights), engine}; } - /// Overload that returns a view closure to compose with other views. + /// Sample from any range. /** - * \param n The number of samples to take. + * If the input range is a particle range, it will extract the weights and treat it as a weighted range. + * The new particles will all have a weight equal to 1, since, after resampling, the probability will be + * represented by the number of particles rather than their individual weight. + * + * If the input range is not a particle range, it will assume a uniform distribution. */ - constexpr auto operator()(std::size_t n) const { - return ranges::make_view_closure(ranges::bind_back(low_variance_sample_fn{}, n)); + template + constexpr auto low_variance_sample_from_range(Range&& range, URNG& engine) const { + static_assert(ranges::sized_range); + static_assert(ranges::random_access_range); + if constexpr (beluga::is_particle_range_v) { + return low_variance_sample_from_range(beluga::views::states(range), beluga::views::weights(range), engine) | + ranges::views::transform(beluga::make_from_state>); + } else { + // Generate uniform weights if no weights are provided. + auto uniform_weights = + ranges::views::repeat_n(1.0, ranges::size(range)); // Use 1.0 for double weights. Important. + return low_variance_sample_from_range(std::forward(range), uniform_weights, engine); + } + } + + // /// Sample from random distributions. + // template + // constexpr auto low_variance_sample_from_distribution(Distribution distribution, URNG& engine) const { + // return ranges::views::generate( + // [distribution = std::move(distribution), &engine]() mutable { return distribution(engine); }); + // } +}; + +/// Implementation detail for a sample range adaptor object. +struct low_variance_sample_fn : public low_variance_sample_base_fn { + /// Overload that takes three arguments. + template + constexpr auto operator()(T&& t, U&& u, V& v) const { + static_assert(ranges::range); + static_assert(ranges::range); + return low_variance_sample_from_range(std::forward(t), std::forward(u), v); // Assume V is a URNG + } + + /// Overload that takes two arguments. + template + constexpr auto operator()(T&& t, U&& u) const { + if constexpr (ranges::range && ranges::range) { + auto& engine = ranges::detail::get_random_engine(); + return low_variance_sample_from_range(std::forward(t), std::forward(u), engine); + // } else if constexpr (is_random_distribution_v) { + // static_assert(std::is_lvalue_reference_v); // Assume U is a URNG + // return low_variance_sample_from_distribution(std::forward(t), u); + } else { + static_assert(ranges::range); + static_assert(std::is_lvalue_reference_v); // Assume U is a URNG + return low_variance_sample_from_range(std::forward(t), u); + } + } + + /// Overload that takes one argument. + template + constexpr auto operator()(T&& t) const { + if constexpr (ranges::range) { + auto& engine = ranges::detail::get_random_engine(); + return low_variance_sample_from_range(std::forward(t), engine); + // } else if constexpr (is_random_distribution_v) { + // auto& engine = ranges::detail::get_random_engine(); + // return low_variance_sample_from_distribution(std::forward(t), engine); + } else { + static_assert(std::is_lvalue_reference_v); // Assume T is a URNG + return ranges::make_view_closure(ranges::bind_back(low_variance_sample_fn{}, std::ref(t))); + } + } + + /// Overload that unwraps the engine reference from a view closure. + template + constexpr auto operator()(Range&& range, std::reference_wrapper engine) const { + static_assert(ranges::range); + return low_variance_sample_from_range(std::forward(range), engine.get()); } }; } // namespace detail /// [Range adaptor object](https://en.cppreference.com/w/cpp/named_req/RangeAdaptorObject) that -/// will perform low-variance sampling on the input range using the provided algorithm. +/// will randomly sample with replacement from an input range. /** - * The input range should be a set of particles with associated weights. - * This adaptor will return a new range of resampled particles using the low-variance sampling - * algorithm described. + * Unlike `std::views::sample`, this does not require a size parameter and the samples will be taken + * from the population **with replacement**, making the sample values independent. + * To use this, the input range must model the + * [random_access_range](https://en.cppreference.com/w/cpp/ranges/random_access_range) + * and [sized_range](https://en.cppreference.com/w/cpp/ranges/sized_range) concepts. + * + * This view implements multinomial resampling for a given range of particles. + * The core idea is to draw random indices / iterators to the input particle range + * from a [multinomial distribution](https://en.wikipedia.org/wiki/Multinomial_distribution) + * parameterized after particle weights (and assumed uniform for non-weighted particle ranges). + * + * This view can also be used to convert any random distribution (a callable that takes a URNG as an + * input argument) into an infinite view that generates values from that distribution. + * + * This view is not cheap to copy, so care must be taken when moving it around. + * Range-v3 does not support move-only views at the time of this implementation. */ -inline constexpr detail::low_variance_sample_fn low_variance_sample; - -} // namespace views +inline constexpr ranges::views::view_closure low_variance_sample; -} // namespace beluga +} // namespace beluga::views #endif From 015f6130647bd13db73004d41bf0610eb88847f5 Mon Sep 17 00:00:00 2001 From: ADEGA Date: Mon, 12 May 2025 15:45:00 -0500 Subject: [PATCH 3/8] [Refactor] Adding algorithm description Signed-off-by: ADEGA --- .../beluga/views/low_variance_sample.hpp | 87 ++++---- beluga/test/beluga/CMakeLists.txt | 1 + .../beluga/views/test_low_variance_sample.cpp | 187 ++++++++++++++++++ 3 files changed, 221 insertions(+), 54 deletions(-) create mode 100644 beluga/test/beluga/views/test_low_variance_sample.cpp diff --git a/beluga/include/beluga/views/low_variance_sample.hpp b/beluga/include/beluga/views/low_variance_sample.hpp index 99e8105ba7..5df44df53f 100644 --- a/beluga/include/beluga/views/low_variance_sample.hpp +++ b/beluga/include/beluga/views/low_variance_sample.hpp @@ -1,4 +1,4 @@ -// Copyright 2024 Ekumen, Inc. +// Copyright 2025 Ekumen, Inc. // // Licensed under the Apache License, Version 2.0 (the "License"); // you may not use this file except in compliance with the License. @@ -20,6 +20,8 @@ #include #include #include +#include +#include #include #include @@ -33,7 +35,16 @@ namespace beluga::views { namespace detail { -/// Implementation of the sample view. +/// Implementation of the low variance sample view. +/// This algorithm computes a single random number r in the range [0, 1/M) and then selects samples from +/// according to this number but still with a probability proportional to the sample weight. This is +/// done by drawing a random number r in the interval [0;M^-1]. Where M is the number of samples in the range. +/// It selects the particles by repeadly adding the fixed amount M^(-1) to r and by choosing the particle that +/// corresponds to to the resulting number. +/// Where M is the number of samples in the range. +/// [Based on] Sebastian Thrun, Wolfram Burgard, and Dieter Fox. 2005. Probabilistic Robotics (Intelligent Robotics +/// and Autonomous Agents). The MIT Press. + /** * \tparam Range A [random access](https://en.cppreference.com/w/cpp/ranges/random_access_range) and * [sized](https://en.cppreference.com/w/cpp/ranges/sized_range) range. @@ -77,41 +88,42 @@ struct low_variance_sample_view /// Construct a cursor from the parent view elements. constexpr explicit cursor(low_variance_sample_view* view) : view_(view), - range_begin_{ranges::begin(view_->range_)}, - weights_begin_{ranges::begin(view_->weights_)}, - M_{static_cast(ranges::size(view_->range_))}, + range_begin_{ranges::begin(view_->range_)}, // Begin iterator of the range + it_{range_begin_}, + weights_begin_{ranges::begin(view_->weights_)}, // Begin iterator of the weights + M_{static_cast(ranges::size(view_->range_))}, // Number of particles to be sampled r_{std::uniform_real_distribution(0.0, 1.0 / static_cast(M_))( - *view_->engine_)}, // Generate r only once - m_{0}, - i_{0}, - c_{*weights_begin_} {} + *view_->engine_)}, // random number in [0, 1/M) + m_{0}, // Current sample index + i_{0}, // Current index in the range. + c_{*weights_begin_} // Cumulative weight of the first particle + {} /// Access the current iterator. - [[nodiscard]] constexpr decltype(auto) read() const noexcept(noexcept(*range_begin_)) { return *it_; } + [[nodiscard]] constexpr decltype(auto) read() const noexcept(noexcept(*this->it_)) { return *it_; } /// Position the current iterator. constexpr void next() { ++m_; + // A number U in [0, 1] that points to exactly one particle in the range. + // Where the particle i satisfies with i=argmin_j \sum_1^j w^m >= U. const double U = r_ + (static_cast(m_) - 1.0) / static_cast(M_); - while (U > c_) { + while (i_ < M_ - 1 && U > c_) { ++i_; c_ += static_cast(*std::next(weights_begin_, i_)); // weights_begin_[i_] } it_ = std::next(range_begin_, i_); // range_begin_ + i_ } - /// Returns true if the cursor is at the end of the range. - constexpr bool equal(const cursor& other) const noexcept { return m_ == other.M_; } - private: low_variance_sample_view* view_; ranges::iterator_t range_begin_; ranges::iterator_t it_; ranges::iterator_t weights_begin_; - long long M_; + uint64_t M_; double r_; - long long m_; - long long i_; + uint64_t m_; + uint64_t i_; double c_; }; @@ -155,18 +167,10 @@ struct low_variance_sample_base_fn { ranges::views::transform(beluga::make_from_state>); } else { // Generate uniform weights if no weights are provided. - auto uniform_weights = - ranges::views::repeat_n(1.0, ranges::size(range)); // Use 1.0 for double weights. Important. + auto uniform_weights = ranges::views::repeat(1.0) | ranges::views::take(ranges::size(range)); return low_variance_sample_from_range(std::forward(range), uniform_weights, engine); } } - - // /// Sample from random distributions. - // template - // constexpr auto low_variance_sample_from_distribution(Distribution distribution, URNG& engine) const { - // return ranges::views::generate( - // [distribution = std::move(distribution), &engine]() mutable { return distribution(engine); }); - // } }; /// Implementation detail for a sample range adaptor object. @@ -185,12 +189,9 @@ struct low_variance_sample_fn : public low_variance_sample_base_fn { if constexpr (ranges::range && ranges::range) { auto& engine = ranges::detail::get_random_engine(); return low_variance_sample_from_range(std::forward(t), std::forward(u), engine); - // } else if constexpr (is_random_distribution_v) { - // static_assert(std::is_lvalue_reference_v); // Assume U is a URNG - // return low_variance_sample_from_distribution(std::forward(t), u); } else { static_assert(ranges::range); - static_assert(std::is_lvalue_reference_v); // Assume U is a URNG + static_assert(std::is_lvalue_reference_v); return low_variance_sample_from_range(std::forward(t), u); } } @@ -201,11 +202,8 @@ struct low_variance_sample_fn : public low_variance_sample_base_fn { if constexpr (ranges::range) { auto& engine = ranges::detail::get_random_engine(); return low_variance_sample_from_range(std::forward(t), engine); - // } else if constexpr (is_random_distribution_v) { - // auto& engine = ranges::detail::get_random_engine(); - // return low_variance_sample_from_distribution(std::forward(t), engine); } else { - static_assert(std::is_lvalue_reference_v); // Assume T is a URNG + static_assert(std::is_lvalue_reference_v); return ranges::make_view_closure(ranges::bind_back(low_variance_sample_fn{}, std::ref(t))); } } @@ -219,27 +217,8 @@ struct low_variance_sample_fn : public low_variance_sample_base_fn { }; } // namespace detail +/// \brief A view adaptor that samples elements from a range with replacement using low variance sampling. -/// [Range adaptor object](https://en.cppreference.com/w/cpp/named_req/RangeAdaptorObject) that -/// will randomly sample with replacement from an input range. -/** - * Unlike `std::views::sample`, this does not require a size parameter and the samples will be taken - * from the population **with replacement**, making the sample values independent. - * To use this, the input range must model the - * [random_access_range](https://en.cppreference.com/w/cpp/ranges/random_access_range) - * and [sized_range](https://en.cppreference.com/w/cpp/ranges/sized_range) concepts. - * - * This view implements multinomial resampling for a given range of particles. - * The core idea is to draw random indices / iterators to the input particle range - * from a [multinomial distribution](https://en.wikipedia.org/wiki/Multinomial_distribution) - * parameterized after particle weights (and assumed uniform for non-weighted particle ranges). - * - * This view can also be used to convert any random distribution (a callable that takes a URNG as an - * input argument) into an infinite view that generates values from that distribution. - * - * This view is not cheap to copy, so care must be taken when moving it around. - * Range-v3 does not support move-only views at the time of this implementation. - */ inline constexpr ranges::views::view_closure low_variance_sample; } // namespace beluga::views diff --git a/beluga/test/beluga/CMakeLists.txt b/beluga/test/beluga/CMakeLists.txt index b7a2103bd6..dfa93c9ba8 100644 --- a/beluga/test/beluga/CMakeLists.txt +++ b/beluga/test/beluga/CMakeLists.txt @@ -67,6 +67,7 @@ add_executable( type_traits/test_tuple_traits.cpp utility/test_forward_like.cpp utility/test_indexing_iterator.cpp + views/test_low_variance_sample.cpp views/test_random_intersperse.cpp views/test_sample.cpp views/test_take_evenly.cpp diff --git a/beluga/test/beluga/views/test_low_variance_sample.cpp b/beluga/test/beluga/views/test_low_variance_sample.cpp new file mode 100644 index 0000000000..037581bb92 --- /dev/null +++ b/beluga/test/beluga/views/test_low_variance_sample.cpp @@ -0,0 +1,187 @@ +// Copyright 2025 Ekumen, Inc. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "beluga/containers/tuple_vector.hpp" +#include "beluga/primitives.hpp" +#include "beluga/views/low_variance_sample.hpp" +#include "beluga/views/particles.hpp" + +namespace { + +TEST(LowVarianceSampleView, FromEmptyRange) { + auto input = std::vector{}; + ASSERT_DEBUG_DEATH(beluga::views::low_variance_sample(input), "Assertion"); +} + +TEST(LowVarianceSampleView, ConceptChecksFromContiguousRange) { + auto input = std::array{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; + auto output = beluga::views::low_variance_sample(input); + + static_assert(ranges::common_range); + static_assert(!ranges::common_range); + + static_assert(!ranges::viewable_range); + static_assert(ranges::viewable_range); + + static_assert(ranges::forward_range); + static_assert(!ranges::forward_range); + + static_assert(ranges::sized_range); + static_assert(!ranges::sized_range); + + static_assert(ranges::bidirectional_range); + static_assert(!ranges::bidirectional_range); + + static_assert(ranges::random_access_range); + static_assert(!ranges::random_access_range); + + static_assert(ranges::contiguous_range); + static_assert(!ranges::contiguous_range); + + static_assert(ranges::range); + static_assert(ranges::semiregular); + static_assert(ranges::enable_view); +} + +TEST(LowVarianceSampleView, UniformDistributionSingleElement) { + auto input = std::array{5}; + auto output = input | beluga::views::low_variance_sample | ranges::views::take_exactly(20); + ASSERT_EQ(ranges::count(output, 5), 20); +} + +TEST(LowVarianceSampleView, DiscreteDistributionSingleElement) { + auto input = std::array{5}; + auto weights = std::array{1.0}; + auto output = beluga::views::low_variance_sample(input, weights) | ranges::views::take_exactly(20); + ASSERT_EQ(ranges::count(output, 5), 20); +} + +TEST(LowVarianceSampleView, DiscreteDistributionSingleElementFromParticleRange) { + auto input = std::array{std::make_tuple(5, beluga::Weight(5.0))}; + auto output = input | beluga::views::low_variance_sample | ranges::views::take_exactly(20) | ranges::to; + ASSERT_EQ(ranges::count(output | beluga::views::states, 5), 20); + ASSERT_EQ(ranges::count(output | beluga::views::weights, beluga::Weight(1.0)), 20); +} + +TEST(LowVarianceSampleView, DoubleDereference) { + auto engine = std::mt19937{std::random_device()()}; + auto input = std::array{10, 42, 39, 20, 50}; + auto output = beluga::views::low_variance_sample(input, engine); + auto it = ranges::begin(output); + ++it; + auto value = *it; + ASSERT_EQ(value, *it); + ASSERT_EQ(value, *it); +} + +TEST(LowVarianceSampleView, NonBorrowedRange) { + auto input = std::array{42}; + const auto create_view = [&]() { return input | beluga::views::low_variance_sample; }; + auto it = ranges::find(create_view(), 42); + static_assert(std::is_same_v); +} + +TEST(LowVarianceSampleView, EngineArgument) { + auto engine = std::mt19937{std::random_device()()}; + auto input = std::array{5}; + auto weights = std::array{1.0}; + auto particles = std::array{std::make_tuple(5, beluga::Weight(1.0))}; + [[maybe_unused]] auto view1 = beluga::views::low_variance_sample(engine); + [[maybe_unused]] auto view2 = input | view1; + [[maybe_unused]] auto view3 = particles | view1; + [[maybe_unused]] auto view4 = beluga::views::low_variance_sample(input, engine); + [[maybe_unused]] auto view5 = beluga::views::low_variance_sample(input, weights, engine); + [[maybe_unused]] auto view6 = beluga::views::low_variance_sample(particles, engine); +} + +// Additional tests specific to low variance sampling properties + +TEST(LowVarianceSampleView, DeterministicWithSameSeed) { + const auto size = 100; + + const auto input = beluga::TupleVector>{ + std::make_tuple(1, beluga::Weight(0.3)), // + std::make_tuple(2, beluga::Weight(0.7)) // + }; + + // First run + auto engine1 = std::mt19937{12345}; + auto output1 = + input | beluga::views::low_variance_sample(engine1) | ranges::views::take_exactly(size) | ranges::to; + + // Second run with same seed + auto engine2 = std::mt19937{12345}; + auto output2 = + input | beluga::views::low_variance_sample(engine2) | ranges::views::take_exactly(size) | ranges::to; + + // Results should be identical + ASSERT_EQ(output1.size(), output2.size()); + for (std::size_t i = 0; i < output1.size(); ++i) { + ASSERT_EQ(beluga::views::states(output1)[i], beluga::views::states(output2)[i]); + } +} + +TEST(LowVarianceSampleView, VerySmallWeights) { + auto input = std::array{1, 2, 3}; + auto weights = std::array{1e-10, 1e-8, 1.0}; + auto output = + beluga::views::low_variance_sample(input, weights) | ranges::views::take_exactly(100) | ranges::to; + + // Element 3 should dominate due to much larger weight + auto count_3 = ranges::count(output, 3); + ASSERT_GT(count_3, 95); // Should be almost all element 3 +} + +TEST(LowVarianceSampleView, LargeRangePerformance) { + // Test with larger input to verify performance characteristics + std::vector large_input(1000); + std::iota(large_input.begin(), large_input.end(), 1); + + std::vector uniform_weights(1000, 1.0 / 1000.0); + + auto engine = std::mt19937{42}; + auto output = beluga::views::low_variance_sample(large_input, uniform_weights, engine) | + ranges::views::take_exactly(5000) | ranges::to; + + // Should complete without issues and produce correct number of samples + ASSERT_EQ(output.size(), 5000); + + // All samples should be from the original range + for (const auto& sample : output) { + ASSERT_GE(sample, 1); + ASSERT_LE(sample, 1000); + } +} + +} // namespace From 61e220b8f37553db03966508dfda3dc082ab28a5 Mon Sep 17 00:00:00 2001 From: ADEGA Date: Mon, 23 Jun 2025 16:43:06 -0500 Subject: [PATCH 4/8] Addressing comments and adding microbenchmark Signed-off-by: ADEGA --- .../beluga/views/low_variance_sample.hpp | 26 +-- .../beluga/views/test_low_variance_sample.cpp | 93 ++++++++- beluga/test/benchmark/CMakeLists.txt | 1 + .../benchmark_low_variance_sample.cpp | 193 ++++++++++++++++++ 4 files changed, 300 insertions(+), 13 deletions(-) create mode 100644 beluga/test/benchmark/benchmark_low_variance_sample.cpp diff --git a/beluga/include/beluga/views/low_variance_sample.hpp b/beluga/include/beluga/views/low_variance_sample.hpp index 5df44df53f..52e9cc48ee 100644 --- a/beluga/include/beluga/views/low_variance_sample.hpp +++ b/beluga/include/beluga/views/low_variance_sample.hpp @@ -36,14 +36,15 @@ namespace beluga::views { namespace detail { /// Implementation of the low variance sample view. -/// This algorithm computes a single random number r in the range [0, 1/M) and then selects samples from -/// according to this number but still with a probability proportional to the sample weight. This is -/// done by drawing a random number r in the interval [0;M^-1]. Where M is the number of samples in the range. -/// It selects the particles by repeadly adding the fixed amount M^(-1) to r and by choosing the particle that -/// corresponds to to the resulting number. -/// Where M is the number of samples in the range. -/// [Based on] Sebastian Thrun, Wolfram Burgard, and Dieter Fox. 2005. Probabilistic Robotics (Intelligent Robotics -/// and Autonomous Agents). The MIT Press. +/// This algorithm is designed to sample elements from a range with replacement, ensuring that the probability of +/// selecting each element is proportional to its weight. It works by computing a single random number r is generated in +/// the range [0, 1/M), where M is the total number of elements in the range. Then, starting from this random number r, +/// the algorithm repeatedly adds a fixed step size of 1/M to r. This step size ensures that the sampling process +/// progresses evenly across the range. After that, for each resulting value, the algorithm selects the element whose +/// cumulative weight corresponds to the current value of r. This ensures that elements with higher weights are more +/// likely to be selected. Finally, the process continues until the desired number of samples is obtained. +/// [Based on] Sebastian Thrun, Wolfram Burgard, and Dieter Fox. 2005. Probabilistic Robotics (Intelligent Robotics and +/// Autonomous Agents). The MIT Press. /** * \tparam Range A [random access](https://en.cppreference.com/w/cpp/ranges/random_access_range) and @@ -104,15 +105,16 @@ struct low_variance_sample_view /// Position the current iterator. constexpr void next() { - ++m_; // A number U in [0, 1] that points to exactly one particle in the range. // Where the particle i satisfies with i=argmin_j \sum_1^j w^m >= U. const double U = r_ + (static_cast(m_) - 1.0) / static_cast(M_); while (i_ < M_ - 1 && U > c_) { ++i_; - c_ += static_cast(*std::next(weights_begin_, i_)); // weights_begin_[i_] + ++weights_begin_; + c_ += *weights_begin_; } - it_ = std::next(range_begin_, i_); // range_begin_ + i_ + ++m_; + it_ = std::next(range_begin_, i_); } private: @@ -124,7 +126,7 @@ struct low_variance_sample_view double r_; uint64_t m_; uint64_t i_; - double c_; + ranges::range_value_t c_; }; /// Return the cursor for the begin iterator. diff --git a/beluga/test/beluga/views/test_low_variance_sample.cpp b/beluga/test/beluga/views/test_low_variance_sample.cpp index 037581bb92..df37404746 100644 --- a/beluga/test/beluga/views/test_low_variance_sample.cpp +++ b/beluga/test/beluga/views/test_low_variance_sample.cpp @@ -168,7 +168,7 @@ TEST(LowVarianceSampleView, LargeRangePerformance) { std::vector large_input(1000); std::iota(large_input.begin(), large_input.end(), 1); - std::vector uniform_weights(1000, 1.0 / 1000.0); + std::vector uniform_weights(large_input.size(), 1.0 / static_cast(large_input.size())); auto engine = std::mt19937{42}; auto output = beluga::views::low_variance_sample(large_input, uniform_weights, engine) | @@ -184,4 +184,95 @@ TEST(LowVarianceSampleView, LargeRangePerformance) { } } +// Multinomial sampling for comparison - simple implementation +std::vector multinomial_sample( + const std::vector& input, + const std::vector& weights, + int num_samples, + std::mt19937& engine) { + std::discrete_distribution<> dist(weights.begin(), weights.end()); + std::vector result; + result.reserve(num_samples); + + for (int i = 0; i < num_samples; ++i) { + result.push_back(input[dist(engine)]); + } + + return result; +} + +TEST(LowVarianceSampleView, LowerVarianceThanMultinomialSampling) { + const auto num_samples = 1000; + const auto num_trials = 200; + + auto input = std::vector{1, 2, 3, 4}; + auto weights = std::vector{0.1, 0.25, 0.35, 0.3}; + + std::vector lv_variances, multinomial_variances; + + for (int trial = 0; trial < num_trials; ++trial) { + // Low variance sampling + auto engine_lv = std::mt19937{static_cast(trial + 5000)}; + auto lv_output = beluga::views::low_variance_sample(input, weights, engine_lv) | + ranges::views::take_exactly(num_samples) | ranges::to; + + // Multinomial sampling + auto engine_mult = std::mt19937{static_cast(trial + 5000)}; // Same seed + auto mult_output = multinomial_sample(input, weights, num_samples, engine_mult); + + auto lv_count_1 = ranges::count(lv_output, 1); + auto mult_count_1 = ranges::count(mult_output, 1); + + double lv_ratio = static_cast(lv_count_1) / num_samples; + double mult_ratio = static_cast(mult_count_1) / num_samples; + + // Store squared deviation from expected (0.1) + lv_variances.push_back(std::pow(lv_ratio - 0.1, 2)); + multinomial_variances.push_back(std::pow(mult_ratio - 0.1, 2)); + } + + // Calculate mean squared errors (approximation of variance) + double lv_mse = std::accumulate(lv_variances.begin(), lv_variances.end(), 0.0) / num_trials; + double mult_mse = std::accumulate(multinomial_variances.begin(), multinomial_variances.end(), 0.0) / num_trials; + + // Low variance sampling should have lower variance (MSE) + ASSERT_LT(lv_mse, mult_mse); + + ASSERT_LT(lv_mse, mult_mse * 0.8); + + std::vector lv_total_var, mult_total_var; + + for (int trial = 0; trial < std::min(50, num_trials); ++trial) { // Subset for performance + auto engine_lv = std::mt19937{static_cast(trial + 6000)}; + auto lv_output = beluga::views::low_variance_sample(input, weights, engine_lv) | + ranges::views::take_exactly(num_samples) | ranges::to; + + auto engine_mult = std::mt19937{static_cast(trial + 6000)}; + auto mult_output = multinomial_sample(input, weights, num_samples, engine_mult); + + // Calculate total variance across all elements + double lv_var_sum = 0.0, mult_var_sum = 0.0; + for (size_t i = 0; i < input.size(); ++i) { + auto lv_count = ranges::count(lv_output, input[i]); + auto mult_count = ranges::count(mult_output, input[i]); + + double lv_ratio = static_cast(lv_count) / num_samples; + double mult_ratio = static_cast(mult_count) / num_samples; + + lv_var_sum += std::pow(lv_ratio - weights[i], 2); + mult_var_sum += std::pow(mult_ratio - weights[i], 2); + } + + lv_total_var.push_back(lv_var_sum); + mult_total_var.push_back(mult_var_sum); + } + + double lv_avg_total_var = + std::accumulate(lv_total_var.begin(), lv_total_var.end(), 0.0) / static_cast(lv_total_var.size()); + double mult_avg_total_var = + std::accumulate(mult_total_var.begin(), mult_total_var.end(), 0.0) / static_cast(mult_total_var.size()); + + ASSERT_LT(lv_avg_total_var, mult_avg_total_var); +} + } // namespace diff --git a/beluga/test/benchmark/CMakeLists.txt b/beluga/test/benchmark/CMakeLists.txt index 209c35065b..9a59634d3a 100644 --- a/beluga/test/benchmark/CMakeLists.txt +++ b/beluga/test/benchmark/CMakeLists.txt @@ -18,6 +18,7 @@ option(BELUGA_RUN_PERFORMANCE_TESTS add_executable( benchmark_beluga benchmark_likelihood_field_model.cpp + benchmark_low_variance_sample.cpp benchmark_main.cpp benchmark_raycasting.cpp benchmark_spatial_hash.cpp diff --git a/beluga/test/benchmark/benchmark_low_variance_sample.cpp b/beluga/test/benchmark/benchmark_low_variance_sample.cpp new file mode 100644 index 0000000000..29065295a9 --- /dev/null +++ b/beluga/test/benchmark/benchmark_low_variance_sample.cpp @@ -0,0 +1,193 @@ +// Copyright 2025 Ekumen, Inc. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "beluga/containers/tuple_vector.hpp" +#include "beluga/primitives.hpp" +#include "beluga/views/low_variance_sample.hpp" +#include "beluga/views/sample.hpp" + +namespace { + +struct State { + double x = 0.; + double y = 0.; + double theta = 0.; +}; + +using Container = beluga::TupleVector>; +using Particle = typename Container::value_type; + +// Helper function to create a container with uniform weights +Container create_uniform_container(std::size_t size) { + auto container = Container{size}; + for (auto&& [state, weight] : container) { + weight = 1.0; + } + return container; +} + +// Helper function to create a container with varied weights +Container create_weighted_container(std::size_t size) { + auto container = Container{size}; + std::mt19937 gen(42); // Fixed seed for reproducibility + std::uniform_real_distribution dist(0.1, 2.0); + + for (auto&& [state, weight] : container) { + weight = dist(gen); + } + return container; +} + +// Helper function to create a container with exponentially distributed weights +Container create_exponential_container(std::size_t size) { + auto container = Container{size}; + std::mt19937 gen(42); // Fixed seed for reproducibility + std::exponential_distribution dist(1.0); + + for (auto&& [state, weight] : container) { + weight = dist(gen) + 0.01; // Add small offset to avoid zero weights + } + return container; +} + +} // namespace + +// Benchmark low variance sampling with uniform weights +void BM_LowVarianceSample_Uniform(benchmark::State& state) { + const auto particle_count = state.range(0); + state.SetComplexityN(particle_count); + const auto container_size = static_cast(particle_count); + + auto container = create_uniform_container(container_size); + auto new_container = Container{container_size}; + + for (auto _ : state) { + auto samples = container | // + beluga::views::low_variance_sample | // + ranges::views::take_exactly(container_size); + auto first = ranges::begin(new_container); + auto last = ranges::copy(samples, first).out; + auto result = ranges::make_subrange(first, last); + benchmark::DoNotOptimize(result); + } +} + +BENCHMARK(BM_LowVarianceSample_Uniform)->RangeMultiplier(2)->Range(128, 1'000'000)->Complexity(); + +// Benchmark low variance sampling with varied weights +void BM_LowVarianceSample_Weighted(benchmark::State& state) { + const auto particle_count = state.range(0); + state.SetComplexityN(particle_count); + const auto container_size = static_cast(particle_count); + + auto container = create_weighted_container(container_size); + auto new_container = Container{container_size}; + + for (auto _ : state) { + auto samples = container | // + beluga::views::low_variance_sample | // + ranges::views::take_exactly(container_size); + auto first = ranges::begin(new_container); + auto last = ranges::copy(samples, first).out; + auto result = ranges::make_subrange(first, last); + benchmark::DoNotOptimize(result); + } +} + +BENCHMARK(BM_LowVarianceSample_Weighted)->RangeMultiplier(2)->Range(128, 1'000'000)->Complexity(); + +// Benchmark low variance sampling with exponentially distributed weights +void BM_LowVarianceSample_Exponential(benchmark::State& state) { + const auto particle_count = state.range(0); + state.SetComplexityN(particle_count); + const auto container_size = static_cast(particle_count); + + auto container = create_exponential_container(container_size); + auto new_container = Container{container_size}; + + for (auto _ : state) { + auto samples = container | // + beluga::views::low_variance_sample | // + ranges::views::take_exactly(container_size); + auto first = ranges::begin(new_container); + auto last = ranges::copy(samples, first).out; + auto result = ranges::make_subrange(first, last); + benchmark::DoNotOptimize(result); + } +} + +BENCHMARK(BM_LowVarianceSample_Exponential)->RangeMultiplier(2)->Range(128, 1'000'000)->Complexity(); + +// Benchmark sampling performance with different sample sizes (fixed population) +void BM_LowVarianceSample_VariableSampleSize(benchmark::State& state) { + const auto sample_count = state.range(0); + state.SetComplexityN(sample_count); + const auto container_size = 10000; // Fixed large population + const auto sample_size = static_cast(sample_count); + + auto container = create_weighted_container(container_size); + auto new_container = Container{sample_size}; + + for (auto _ : state) { + auto samples = container | // + beluga::views::low_variance_sample | // + ranges::views::take_exactly(sample_size); + auto first = ranges::begin(new_container); + auto last = ranges::copy(samples, first).out; + auto result = ranges::make_subrange(first, last); + benchmark::DoNotOptimize(result); + state.counters["SampleSize"] = static_cast(sample_size); + state.counters["PopulationSize"] = static_cast(container_size); + } +} + +BENCHMARK(BM_LowVarianceSample_VariableSampleSize)->RangeMultiplier(2)->Range(64, 8192)->Complexity(); + +// Benchmark memory access patterns by testing small vs large sample sizes +void BM_LowVarianceSample_MemoryPattern(benchmark::State& state) { + const auto particle_count = state.range(0); + state.SetComplexityN(particle_count); + const auto container_size = static_cast(particle_count); + const auto sample_size = std::max(static_cast(1), container_size / 10); // Sample 10% of particles + + auto container = create_weighted_container(container_size); + auto new_container = Container{sample_size}; + + for (auto _ : state) { + auto samples = container | // + beluga::views::low_variance_sample | // + ranges::views::take_exactly(sample_size); + auto first = ranges::begin(new_container); + auto last = ranges::copy(samples, first).out; + auto result = ranges::make_subrange(first, last); + benchmark::DoNotOptimize(result); + state.counters["SampleRatio"] = static_cast(sample_size) / static_cast(container_size); + } +} + +BENCHMARK(BM_LowVarianceSample_MemoryPattern)->RangeMultiplier(2)->Range(1000, 100'000)->Complexity(); From 0ea1a57f333ede97b686cf14593e13c72d9f0286 Mon Sep 17 00:00:00 2001 From: ADEGA Date: Mon, 23 Jun 2025 16:51:42 -0500 Subject: [PATCH 5/8] Adding information about low-variance-sample to readme Signed-off-by: ADEGA --- beluga/docs/_doxygen/beluga-main.md | 1 + beluga/docs/index.md | 1 + 2 files changed, 2 insertions(+) diff --git a/beluga/docs/_doxygen/beluga-main.md b/beluga/docs/_doxygen/beluga-main.md index fb5de4c6f0..e8e5ccb1d4 100644 --- a/beluga/docs/_doxygen/beluga-main.md +++ b/beluga/docs/_doxygen/beluga-main.md @@ -41,6 +41,7 @@ They are lazily evaluated range adaptor objects compatible with the Range-v3 lib | [beluga::views::states](@ref views/particles.hpp) | Produces a view of the states of a range of particles | | [beluga::views::take_evenly](@ref views/take_evenly.hpp) | Returns a range consisting of `count` elements evenly spaced over the source range | | [beluga::views::take_while_kld](@ref views/take_while_kld.hpp) | Take elements from a range while the KLD condition is statisfied | +| [beluga::views::low_variance_sample](@ref views/low_variance_sample.hpp) | Implements low variance resampling for a given range of particles | | [beluga::views::weights](@ref views/particles.hpp) | Produces a view of the weights of a range of particles | | [beluga::views::zip](@ref views/zip.hpp) | Given N ranges, return a new range where the Mth element is a tuple of the Mth elements of all N ranges | diff --git a/beluga/docs/index.md b/beluga/docs/index.md index 53f74b1c22..d63ce753c2 100644 --- a/beluga/docs/index.md +++ b/beluga/docs/index.md @@ -23,6 +23,7 @@ The current set of features includes: - Multinomial resampling from a particle range - [Adaptive KLD resampling][fox2001] - [Selective resampling][grisetti2007], on-motion resampling, and interval resampling policies + - Low Variance sampling - Support for sequential and parallel execution policies - Weighted mean and covariance statistics for pose estimation - Sensor models: From 60d5ebc53a4d71cc309febf50016fe315a454918 Mon Sep 17 00:00:00 2001 From: ADEGA Date: Mon, 1 Sep 2025 18:09:55 -0500 Subject: [PATCH 6/8] Addressing comments Signed-off-by: ADEGA --- beluga/docs/index.md | 3 +- .../beluga/views/low_variance_sample.hpp | 61 +++---- .../beluga/views/test_low_variance_sample.cpp | 170 +++++++++--------- .../benchmark_low_variance_sample.cpp | 18 +- 4 files changed, 129 insertions(+), 123 deletions(-) diff --git a/beluga/docs/index.md b/beluga/docs/index.md index d63ce753c2..7fc0c7d1b5 100644 --- a/beluga/docs/index.md +++ b/beluga/docs/index.md @@ -23,7 +23,7 @@ The current set of features includes: - Multinomial resampling from a particle range - [Adaptive KLD resampling][fox2001] - [Selective resampling][grisetti2007], on-motion resampling, and interval resampling policies - - Low Variance sampling + - [Low Variance sampling][thrun2005] - Support for sequential and parallel execution policies - Weighted mean and covariance statistics for pose estimation - Sensor models: @@ -45,3 +45,4 @@ Beluga is built on top of the following open source libraries: [aos_soa]: https://en.wikipedia.org/wiki/AoS_and_SoA [fox2001]: https://dl.acm.org/doi/10.5555/2980539.2980632 [grisetti2007]: https://doi.org/10.1109/TRO.2006.889486 +[thrun2005]: https://books.google.com.ar/books?id=jtSMEAAAQBAJ diff --git a/beluga/include/beluga/views/low_variance_sample.hpp b/beluga/include/beluga/views/low_variance_sample.hpp index 52e9cc48ee..c4094a5ab3 100644 --- a/beluga/include/beluga/views/low_variance_sample.hpp +++ b/beluga/include/beluga/views/low_variance_sample.hpp @@ -37,9 +37,9 @@ namespace detail { /// Implementation of the low variance sample view. /// This algorithm is designed to sample elements from a range with replacement, ensuring that the probability of -/// selecting each element is proportional to its weight. It works by computing a single random number r is generated in -/// the range [0, 1/M), where M is the total number of elements in the range. Then, starting from this random number r, -/// the algorithm repeatedly adds a fixed step size of 1/M to r. This step size ensures that the sampling process +/// selecting each element is proportional to its weight. It works with a single random number r in +/// the [0, 1/M) interval, where M is the total number of elements in the range. Then, starting from this random number +/// r, the algorithm repeatedly adds a fixed step size of 1/M to r. This step size ensures that the sampling process /// progresses evenly across the range. After that, for each resulting value, the algorithm selects the element whose /// cumulative weight corresponds to the current value of r. This ensures that elements with higher weights are more /// likely to be selected. Finally, the process continues until the desired number of samples is obtained. @@ -96,7 +96,6 @@ struct low_variance_sample_view r_{std::uniform_real_distribution(0.0, 1.0 / static_cast(M_))( *view_->engine_)}, // random number in [0, 1/M) m_{0}, // Current sample index - i_{0}, // Current index in the range. c_{*weights_begin_} // Cumulative weight of the first particle {} @@ -107,14 +106,13 @@ struct low_variance_sample_view constexpr void next() { // A number U in [0, 1] that points to exactly one particle in the range. // Where the particle i satisfies with i=argmin_j \sum_1^j w^m >= U. - const double U = r_ + (static_cast(m_) - 1.0) / static_cast(M_); - while (i_ < M_ - 1 && U > c_) { - ++i_; + const double U = r_ + static_cast(m_) / static_cast(M_); + while (it_ != ranges::end(view_->range_) && U > c_) { + ++it_; ++weights_begin_; c_ += *weights_begin_; } ++m_; - it_ = std::next(range_begin_, i_); } private: @@ -125,7 +123,6 @@ struct low_variance_sample_view uint64_t M_; double r_; uint64_t m_; - uint64_t i_; ranges::range_value_t c_; }; @@ -178,43 +175,43 @@ struct low_variance_sample_base_fn { /// Implementation detail for a sample range adaptor object. struct low_variance_sample_fn : public low_variance_sample_base_fn { /// Overload that takes three arguments. - template - constexpr auto operator()(T&& t, U&& u, V& v) const { - static_assert(ranges::range); - static_assert(ranges::range); - return low_variance_sample_from_range(std::forward(t), std::forward(u), v); // Assume V is a URNG + template + constexpr auto operator()(Range&& range, Weights&& weights, RandomEngine& random_engine) const { + static_assert(ranges::range); + static_assert(ranges::range); + return low_variance_sample_from_range(std::forward(range), std::forward(weights), random_engine); } /// Overload that takes two arguments. - template - constexpr auto operator()(T&& t, U&& u) const { - if constexpr (ranges::range && ranges::range) { - auto& engine = ranges::detail::get_random_engine(); - return low_variance_sample_from_range(std::forward(t), std::forward(u), engine); + template + constexpr auto operator()(Range&& range, Weights&& weights) const { + if constexpr (ranges::range && ranges::range) { + auto& random_engine = ranges::detail::get_random_engine(); + return low_variance_sample_from_range(std::forward(range), std::forward(weights), random_engine); } else { - static_assert(ranges::range); - static_assert(std::is_lvalue_reference_v); - return low_variance_sample_from_range(std::forward(t), u); + static_assert(ranges::range); + static_assert(std::is_lvalue_reference_v); + return low_variance_sample_from_range(std::forward(range), weights); } } /// Overload that takes one argument. - template - constexpr auto operator()(T&& t) const { - if constexpr (ranges::range) { - auto& engine = ranges::detail::get_random_engine(); - return low_variance_sample_from_range(std::forward(t), engine); + template + constexpr auto operator()(Range&& range) const { + if constexpr (ranges::range) { + auto& random_engine = ranges::detail::get_random_engine(); + return low_variance_sample_from_range(std::forward(range), random_engine); } else { - static_assert(std::is_lvalue_reference_v); - return ranges::make_view_closure(ranges::bind_back(low_variance_sample_fn{}, std::ref(t))); + static_assert(std::is_lvalue_reference_v); + return ranges::make_view_closure(ranges::bind_back(low_variance_sample_fn{}, std::ref(range))); } } /// Overload that unwraps the engine reference from a view closure. - template - constexpr auto operator()(Range&& range, std::reference_wrapper engine) const { + template + constexpr auto operator()(Range&& range, std::reference_wrapper random_engine) const { static_assert(ranges::range); - return low_variance_sample_from_range(std::forward(range), engine.get()); + return low_variance_sample_from_range(std::forward(range), random_engine.get()); } }; diff --git a/beluga/test/beluga/views/test_low_variance_sample.cpp b/beluga/test/beluga/views/test_low_variance_sample.cpp index df37404746..9b640cf844 100644 --- a/beluga/test/beluga/views/test_low_variance_sample.cpp +++ b/beluga/test/beluga/views/test_low_variance_sample.cpp @@ -36,6 +36,7 @@ #include "beluga/primitives.hpp" #include "beluga/views/low_variance_sample.hpp" #include "beluga/views/particles.hpp" +#include "beluga/views/sample.hpp" namespace { @@ -80,19 +81,21 @@ TEST(LowVarianceSampleView, UniformDistributionSingleElement) { ASSERT_EQ(ranges::count(output, 5), 20); } -TEST(LowVarianceSampleView, DiscreteDistributionSingleElement) { - auto input = std::array{5}; - auto weights = std::array{1.0}; - auto output = beluga::views::low_variance_sample(input, weights) | ranges::views::take_exactly(20); - ASSERT_EQ(ranges::count(output, 5), 20); -} - -TEST(LowVarianceSampleView, DiscreteDistributionSingleElementFromParticleRange) { - auto input = std::array{std::make_tuple(5, beluga::Weight(5.0))}; - auto output = input | beluga::views::low_variance_sample | ranges::views::take_exactly(20) | ranges::to; - ASSERT_EQ(ranges::count(output | beluga::views::states, 5), 20); - ASSERT_EQ(ranges::count(output | beluga::views::weights, beluga::Weight(1.0)), 20); -} +// TEST(LowVarianceSampleView, DiscreteDistributionSingleElement) { +// auto input = std::array{5}; +// auto weights = std::array{1.0}; +// auto output = beluga::views::low_variance_sample(input, weights) | // +// ranges::views::take_exactly(20) | // +// ranges::to; +// ASSERT_EQ(ranges::count(output, 5), 20); +// } + +// TEST(LowVarianceSampleView, DiscreteDistributionSingleElementFromParticleRange) { +// auto input = std::array{std::make_tuple(5, beluga::Weight(5.0))}; +// auto output = input | beluga::views::low_variance_sample | ranges::views::take_exactly(20) | +// ranges::to; ASSERT_EQ(ranges::count(output | beluga::views::states, 5), 20); +// ASSERT_EQ(ranges::count(output | beluga::views::weights, beluga::Weight(1.0)), 20); +// } TEST(LowVarianceSampleView, DoubleDereference) { auto engine = std::mt19937{std::random_device()()}; @@ -152,27 +155,14 @@ TEST(LowVarianceSampleView, DeterministicWithSameSeed) { } } -TEST(LowVarianceSampleView, VerySmallWeights) { - auto input = std::array{1, 2, 3}; - auto weights = std::array{1e-10, 1e-8, 1.0}; - auto output = - beluga::views::low_variance_sample(input, weights) | ranges::views::take_exactly(100) | ranges::to; - - // Element 3 should dominate due to much larger weight - auto count_3 = ranges::count(output, 3); - ASSERT_GT(count_3, 95); // Should be almost all element 3 -} - TEST(LowVarianceSampleView, LargeRangePerformance) { // Test with larger input to verify performance characteristics std::vector large_input(1000); std::iota(large_input.begin(), large_input.end(), 1); - std::vector uniform_weights(large_input.size(), 1.0 / static_cast(large_input.size())); - auto engine = std::mt19937{42}; - auto output = beluga::views::low_variance_sample(large_input, uniform_weights, engine) | - ranges::views::take_exactly(5000) | ranges::to; + auto output = beluga::views::low_variance_sample(large_input, engine) | ranges::views::take_exactly(5000) | + ranges::to; // Should complete without issues and produce correct number of samples ASSERT_EQ(output.size(), 5000); @@ -184,95 +174,97 @@ TEST(LowVarianceSampleView, LargeRangePerformance) { } } -// Multinomial sampling for comparison - simple implementation -std::vector multinomial_sample( - const std::vector& input, - const std::vector& weights, - int num_samples, - std::mt19937& engine) { - std::discrete_distribution<> dist(weights.begin(), weights.end()); - std::vector result; - result.reserve(num_samples); - - for (int i = 0; i < num_samples; ++i) { - result.push_back(input[dist(engine)]); - } - - return result; -} - -TEST(LowVarianceSampleView, LowerVarianceThanMultinomialSampling) { +TEST(LowVarianceSampleView, LowerVarianceOfSampleMean) { const auto num_samples = 1000; const auto num_trials = 200; auto input = std::vector{1, 2, 3, 4}; auto weights = std::vector{0.1, 0.25, 0.35, 0.3}; - std::vector lv_variances, multinomial_variances; + double true_mean = 0.0; + for (size_t i = 0; i < input.size(); ++i) { + true_mean += static_cast(input[i]) * weights[i]; + } + + std::vector lv_means, sample_means; + lv_means.reserve(num_trials); + sample_means.reserve(num_trials); for (int trial = 0; trial < num_trials; ++trial) { // Low variance sampling auto engine_lv = std::mt19937{static_cast(trial + 5000)}; auto lv_output = beluga::views::low_variance_sample(input, weights, engine_lv) | ranges::views::take_exactly(num_samples) | ranges::to; + const auto lv_sum = std::accumulate(lv_output.begin(), lv_output.end(), 0LL); + lv_means.push_back(static_cast(lv_sum) / num_samples); // Multinomial sampling - auto engine_mult = std::mt19937{static_cast(trial + 5000)}; // Same seed - auto mult_output = multinomial_sample(input, weights, num_samples, engine_mult); - - auto lv_count_1 = ranges::count(lv_output, 1); - auto mult_count_1 = ranges::count(mult_output, 1); + auto engine_sample = std::mt19937{static_cast(trial + 5000)}; // Same seed + auto sample_output = beluga::views::sample(input, weights, engine_sample) | + ranges::views::take_exactly(num_samples) | ranges::to; + const auto sample_sum = std::accumulate(sample_output.begin(), sample_output.end(), 0LL); + sample_means.push_back(static_cast(sample_sum) / num_samples); + } - double lv_ratio = static_cast(lv_count_1) / num_samples; - double mult_ratio = static_cast(mult_count_1) / num_samples; + const auto calculate_variance = [&](const std::vector& means) { + double sum_sq_err = 0.0; + for (double mean : means) { + sum_sq_err += std::pow(mean - true_mean, 2); + } + return sum_sq_err / num_trials; // This is MSE, which is variance for an unbiased estimator + }; - // Store squared deviation from expected (0.1) - lv_variances.push_back(std::pow(lv_ratio - 0.1, 2)); - multinomial_variances.push_back(std::pow(mult_ratio - 0.1, 2)); - } + const double lv_variance = calculate_variance(lv_means); + const double sample_variance = calculate_variance(sample_means); - // Calculate mean squared errors (approximation of variance) - double lv_mse = std::accumulate(lv_variances.begin(), lv_variances.end(), 0.0) / num_trials; - double mult_mse = std::accumulate(multinomial_variances.begin(), multinomial_variances.end(), 0.0) / num_trials; + // Low variance sampling should have lower variance. + ASSERT_LT(lv_variance, sample_variance); +} - // Low variance sampling should have lower variance (MSE) - ASSERT_LT(lv_mse, mult_mse); +TEST(LowVarianceSampleView, LowerSumOfSquaredErrorOfProbabilities) { + const auto num_samples = 1000; + const auto num_trials = 200; - ASSERT_LT(lv_mse, mult_mse * 0.8); + auto input = std::vector{1, 2, 3, 4}; + auto weights = std::vector{0.1, 0.25, 0.35, 0.3}; + const auto limited_trials = std::min(50, num_trials); // Subset for performance - std::vector lv_total_var, mult_total_var; + std::vector lv_total_sse, sample_total_sse; + lv_total_sse.reserve(limited_trials); + sample_total_sse.reserve(limited_trials); - for (int trial = 0; trial < std::min(50, num_trials); ++trial) { // Subset for performance + for (int trial = 0; trial < limited_trials; ++trial) { auto engine_lv = std::mt19937{static_cast(trial + 6000)}; auto lv_output = beluga::views::low_variance_sample(input, weights, engine_lv) | ranges::views::take_exactly(num_samples) | ranges::to; - auto engine_mult = std::mt19937{static_cast(trial + 6000)}; - auto mult_output = multinomial_sample(input, weights, num_samples, engine_mult); - - // Calculate total variance across all elements - double lv_var_sum = 0.0, mult_var_sum = 0.0; - for (size_t i = 0; i < input.size(); ++i) { - auto lv_count = ranges::count(lv_output, input[i]); - auto mult_count = ranges::count(mult_output, input[i]); - - double lv_ratio = static_cast(lv_count) / num_samples; - double mult_ratio = static_cast(mult_count) / num_samples; - - lv_var_sum += std::pow(lv_ratio - weights[i], 2); - mult_var_sum += std::pow(mult_ratio - weights[i], 2); - } - - lv_total_var.push_back(lv_var_sum); - mult_total_var.push_back(mult_var_sum); + auto engine_sample = std::mt19937{static_cast(trial + 6000)}; + auto sample_output = beluga::views::sample(input, weights, engine_sample) | + ranges::views::take_exactly(num_samples) | ranges::to; + + const auto get_sse = [&](const auto& output) { + std::map counts; + for (int value : output) { + counts[value]++; + } + double sse = 0.0; + for (size_t i = 0; i < input.size(); ++i) { + const double ratio = static_cast(counts[input[i]]) / num_samples; + sse += std::pow(ratio - weights[i], 2); + } + return sse; + }; + + lv_total_sse.push_back(get_sse(lv_output)); + sample_total_sse.push_back(get_sse(sample_output)); } - double lv_avg_total_var = - std::accumulate(lv_total_var.begin(), lv_total_var.end(), 0.0) / static_cast(lv_total_var.size()); - double mult_avg_total_var = - std::accumulate(mult_total_var.begin(), mult_total_var.end(), 0.0) / static_cast(mult_total_var.size()); + const auto lv_avg_sse = + std::accumulate(lv_total_sse.begin(), lv_total_sse.end(), 0.0) / static_cast(lv_total_sse.size()); + const auto sample_avg_sse = std::accumulate(sample_total_sse.begin(), sample_total_sse.end(), 0.0) / + static_cast(sample_total_sse.size()); - ASSERT_LT(lv_avg_total_var, mult_avg_total_var); + ASSERT_LT(lv_avg_sse, sample_avg_sse); } } // namespace diff --git a/beluga/test/benchmark/benchmark_low_variance_sample.cpp b/beluga/test/benchmark/benchmark_low_variance_sample.cpp index 29065295a9..176eb39ce5 100644 --- a/beluga/test/benchmark/benchmark_low_variance_sample.cpp +++ b/beluga/test/benchmark/benchmark_low_variance_sample.cpp @@ -57,9 +57,17 @@ Container create_weighted_container(std::size_t size) { std::mt19937 gen(42); // Fixed seed for reproducibility std::uniform_real_distribution dist(0.1, 2.0); + double total_weight = 0.0; for (auto&& [state, weight] : container) { weight = dist(gen); + total_weight += weight; } + + // Normalize weights + for (auto&& [state, weight] : container) { + weight /= total_weight; + } + return container; } @@ -69,9 +77,17 @@ Container create_exponential_container(std::size_t size) { std::mt19937 gen(42); // Fixed seed for reproducibility std::exponential_distribution dist(1.0); + double total_weight = 0.0; for (auto&& [state, weight] : container) { weight = dist(gen) + 0.01; // Add small offset to avoid zero weights + total_weight += weight; } + + // Normalize weights + for (auto&& [state, weight] : container) { + weight /= total_weight; + } + return container; } @@ -173,7 +189,7 @@ void BM_LowVarianceSample_MemoryPattern(benchmark::State& state) { const auto particle_count = state.range(0); state.SetComplexityN(particle_count); const auto container_size = static_cast(particle_count); - const auto sample_size = std::max(static_cast(1), container_size / 10); // Sample 10% of particles + const auto sample_size = std::max(1UL, container_size / 10UL); // Sample 10% of particles auto container = create_weighted_container(container_size); auto new_container = Container{sample_size}; From a41b3d393ffe2264517281234af53e28b7eff7f1 Mon Sep 17 00:00:00 2001 From: "Marco A. Montero" Date: Tue, 21 Apr 2026 14:22:00 -0400 Subject: [PATCH 7/8] Solve CMake lint issue --- beluga/test/beluga/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/beluga/test/beluga/CMakeLists.txt b/beluga/test/beluga/CMakeLists.txt index 8307af947c..50299b3ebd 100644 --- a/beluga/test/beluga/CMakeLists.txt +++ b/beluga/test/beluga/CMakeLists.txt @@ -67,8 +67,8 @@ add_executable( type_traits/test_tuple_traits.cpp utility/test_forward_like.cpp utility/test_indexing_iterator.cpp - views/test_low_variance_sample.cpp views/test_likelihoods.cpp + views/test_low_variance_sample.cpp views/test_random_intersperse.cpp views/test_sample.cpp views/test_take_evenly.cpp From dd95dec62eda28bfe117c2522bb71d24acf4c75a Mon Sep 17 00:00:00 2001 From: "Marco A. Montero" Date: Tue, 21 Apr 2026 14:25:05 -0400 Subject: [PATCH 8/8] Solve include order lint issue --- beluga/include/beluga/views.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/beluga/include/beluga/views.hpp b/beluga/include/beluga/views.hpp index 884b9aa2d0..5fac5bcfb2 100644 --- a/beluga/include/beluga/views.hpp +++ b/beluga/include/beluga/views.hpp @@ -16,8 +16,8 @@ #define BELUGA_VIEWS_HPP #include -#include #include +#include #include #include #include