Implement low-variance sampling algorithm by larodriguez22 · Pull Request #482 · Ekumen-OS/beluga

larodriguez22 · 2025-04-30T03:52:52Z

Proposed changes

This PR adds the strategy of low-variance sampling for the resampling step. Adds a new range adaptor called beluga::views::low_variance_sample. Related to: #48

Type of change

🐛 Bugfix (change which fixes an issue)
🚀 Feature (change which adds functionality)
📚 Documentation (change which fixes or extends documentation)

Checklist

Put an x in the boxes that apply. This is simply a reminder of what we will require before merging your code.

Lint and unit tests (if any) pass locally with my changes
I have added tests that prove my fix is effective or that my feature works
I have added necessary documentation (if appropriate)
All commits have been signed for DCO

Additional comments

How to run it?

Run tests:

colcon build --packages-up-to beluga && ./build/beluga/test/beluga/test_beluga --gtest_filter="*LowVarianceSampleView*"

Run microbenchmarks

colcon build --packages-up-to beluga && ./build/beluga/test/benchmark/benchmark_beluga --benchmark_filter=Low

hidmic

Great first take @larodriguez22! I think adding some unit tests will help a lot.

beluga/include/beluga/views/low_variance_sample.hpp

hidmic

Second pass. Great job @larodriguez22 !

beluga/include/beluga/views/low_variance_sample.hpp

hidmic · 2025-05-13T18:54:31Z

beluga/include/beluga/views/low_variance_sample.hpp

+
+    /// Position the current iterator.
+    constexpr void next() {
+      ++m_;


@larodriguez22 meta: can we use complete names for variables? Or document the algorithm to put u, c, m, M in context?

Yes, I used those names because of the book. But I can change them

For the algorithm itself I don't mind using single letter variables so long as we have a comment block right next to it that explains what those variables mean.

beluga/include/beluga/views/low_variance_sample.hpp

glpuga

Good work @larodriguez22 .

Which of the versions discussed in #48 is targetted here? this looks like the version in Prob Robotics, but that version is known not to work with KLD. Is this meant to be used with fixed resampling?

beluga/include/beluga/views/low_variance_sample.hpp

glpuga · 2025-05-22T13:03:06Z

beluga/include/beluga/views/low_variance_sample.hpp

+          c_{*weights_begin_} {}
+
+    /// Access the current iterator.
+    [[nodiscard]] constexpr decltype(auto) read() const noexcept(noexcept(*range_begin_)) { return *it_; }


noexcept(noexcept(*range_begin_))

C++ needs to burn in hell.

It was like that in the example of SampleView, isn't?

beluga/include/beluga/views/low_variance_sample.hpp

glpuga · 2025-05-22T13:30:18Z

beluga/include/beluga/views/low_variance_sample.hpp

+  template <class T, class U, class V>
+  constexpr auto operator()(T&& t, U&& u, V& v) const {


Let's not use single letter variables and template arguments, and even less if they are unrelated to their role in the code. TUV are just three consecutive letters. Names should be convey meaning. Use whole words, multiple words if needed, we are not billed for the characters we use.

https://google.github.io/styleguide/cppguide.html#General_Naming_Rules

I used those names, because of conventions to the book. But I'll change them. Thanks for the link

This change is pending.

glpuga · 2025-05-22T13:32:09Z

beluga/include/beluga/views/low_variance_sample.hpp

+      static_assert(ranges::range<T>);
+      static_assert(std::is_lvalue_reference_v<U&&>);  // Assume U is a URNG
+      return low_variance_sample_from_range(std::forward<T>(t), u);


I think this needs more explaining.

larodriguez22 · 2025-06-03T14:26:25Z

I created the tests, but @hidmic or @glpuga. Do you have any additional scenarios I could make?

hidmic

Great work @larodriguez22. Second pass.

beluga/include/beluga/views/low_variance_sample.hpp

hidmic · 2025-06-03T21:47:44Z

beluga/include/beluga/views/low_variance_sample.hpp

+    double r_;
+    uint64_t m_;
+    uint64_t i_;
+    double c_;


@larodriguez22 meta: why do we force this to use double types, instead of using the Range value type?

I changed c_ to ranges::range_value_t<Weights>, but I think the other ones must be fixed to be double or int numbers

hidmic · 2025-06-03T21:49:17Z

beluga/include/beluga/views/low_variance_sample.hpp

+
+    /// Position the current iterator.
+    constexpr void next() {
+      ++m_;


For the algorithm itself I don't mind using single letter variables so long as we have a comment block right next to it that explains what those variables mean.

beluga/test/beluga/views/test_low_variance_sample.cpp

hidmic · 2025-06-03T21:52:20Z

log/latest_test

@@ -0,0 +1 @@
+test_2025-06-01_21-06-40


@larodriguez22 let's remove this log directory. I wonder why it was picked up to begin with 🤔

beluga/test/beluga/views/test_low_variance_sample.cpp

beluga/include/beluga/views/low_variance_sample.hpp

hidmic

Third pass. Great work @larodriguez22.

hidmic · 2025-06-24T11:40:01Z

beluga/docs/index.md

  - Multinomial resampling from a particle range
  - [Adaptive KLD resampling][fox2001]
  - [Selective resampling][grisetti2007], on-motion resampling, and interval resampling policies
+  - Low Variance sampling


@larodriguez22 nit: would be nice to add a citation like we do for KLD and selective resampling.

hidmic · 2025-06-24T11:41:17Z

beluga/include/beluga/views/low_variance_sample.hpp

+/// 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,


@larodriguez22 nit:

Suggested change

/// 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,

/// 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,

beluga/include/beluga/views/low_variance_sample.hpp

hidmic · 2025-06-24T11:44:19Z

beluga/include/beluga/views/low_variance_sample.hpp

+        c_ += *weights_begin_;
+      }
+      ++m_;
+      it_ = std::next(range_begin_, i_);


@larodriguez22 do we need i_? Why not increment it_ in the loop, checking it against the end iterator?

You are right, fixed

hidmic · 2025-06-24T11:44:55Z

beluga/include/beluga/views/low_variance_sample.hpp

+  template <class T, class U, class V>
+  constexpr auto operator()(T&& t, U&& u, V& v) const {


This change is pending.

hidmic · 2025-06-24T11:50:36Z

beluga/test/beluga/views/test_low_variance_sample.cpp

+    auto lv_count_1 = ranges::count(lv_output, 1);
+    auto mult_count_1 = ranges::count(mult_output, 1);
+
+    double lv_ratio = static_cast<double>(lv_count_1) / num_samples;
+    double mult_ratio = static_cast<double>(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));


@larodriguez22 this is peculiar. We are computing variances over probabilities. I was expecting variances over distribution statistics (mean, variance, etc.). I wonder if results are equivalent (they should be, right?).

It is not necessarily the same

hidmic · 2025-06-24T11:53:04Z

beluga/test/beluga/views/test_low_variance_sample.cpp

+
+  std::vector<double> lv_total_var, mult_total_var;
+
+  for (int trial = 0; trial < std::min(50, num_trials); ++trial) {  // Subset for performance


@larodriguez22 this should perhaps be a separate test.

hidmic · 2025-06-24T12:08:38Z

beluga/test/benchmark/benchmark_low_variance_sample.cpp

+    auto first = ranges::begin(new_container);
+    auto last = ranges::copy(samples, first).out;
+    auto result = ranges::make_subrange(first, last);


@larodriguez22 why do we do this? Is it to force computation?

Yes, it is to force computation

hidmic · 2025-06-24T12:10:52Z

beluga/test/benchmark/benchmark_low_variance_sample.cpp

+  const auto particle_count = state.range(0);
+  state.SetComplexityN(particle_count);
+  const auto container_size = static_cast<std::size_t>(particle_count);
+  const auto sample_size = std::max(static_cast<std::size_t>(1), container_size / 10);  // Sample 10% of particles


@larodriguez22 nit:

Suggested change

const auto sample_size = std::max(static_cast<std::size_t>(1), container_size / 10); // Sample 10% of particles

const auto sample_size = std::max(1UL, container_size / 10UL); // Sample 10% of particles

hidmic · 2025-06-24T12:15:16Z

beluga/test/benchmark/benchmark_low_variance_sample.cpp

+  for (auto&& [state, weight] : container) {
+    weight = dist(gen);
+  }
+  return container;


@larodriguez22 meta: don't we have to normalize weights? Same elsewhere.

True, fixed

hidmic · 2025-08-26T21:45:42Z

@larodriguez22 may I take this to completion or are you still working on it?

Signed-off-by: ADEGA <la.rodriguez@uniandes.edu.co>

… sample Signed-off-by: ADEGA <la.rodriguez@uniandes.edu.co>

Signed-off-by: ADEGA <la.rodriguez@uniandes.edu.co>

larodriguez22 · 2025-09-01T23:20:21Z

@larodriguez22 may I take this to completion or are you still working on it?

I took longer than expected. I'm sorry, ready for another review

larodriguez22 requested a review from hidmic April 30, 2025 03:52

larodriguez22 force-pushed the larodriguez22/48-low-variance-sampling branch from 4a84380 to 614ed8f Compare April 30, 2025 04:09

hidmic requested changes May 5, 2025

View reviewed changes

larodriguez22 force-pushed the larodriguez22/48-low-variance-sampling branch from bd2804c to 3c63caf Compare May 12, 2025 19:51

hidmic requested changes May 13, 2025

View reviewed changes

glpuga requested changes May 22, 2025

View reviewed changes

hidmic requested review from glpuga and hidmic June 3, 2025 14:17

hidmic requested changes Jun 3, 2025

View reviewed changes

larodriguez22 requested a review from hidmic June 23, 2025 21:43

hidmic reviewed Jun 24, 2025

View reviewed changes

larodriguez22 added 6 commits September 1, 2025 18:13

First start low-variance sampling

b2e8c07

Signed-off-by: ADEGA <la.rodriguez@uniandes.edu.co>

[Refactor] Addressing Comments. Making low_variance_sample similar to…

99341f6

… sample Signed-off-by: ADEGA <la.rodriguez@uniandes.edu.co>

[Refactor] Adding algorithm description

015f613

Signed-off-by: ADEGA <la.rodriguez@uniandes.edu.co>

Addressing comments and adding microbenchmark

61e220b

Signed-off-by: ADEGA <la.rodriguez@uniandes.edu.co>

Adding information about low-variance-sample to readme

0ea1a57

Signed-off-by: ADEGA <la.rodriguez@uniandes.edu.co>

Addressing comments

60d5ebc

Signed-off-by: ADEGA <la.rodriguez@uniandes.edu.co>

larodriguez22 force-pushed the larodriguez22/48-low-variance-sampling branch from 323a77b to 60d5ebc Compare September 1, 2025 23:16

larodriguez22 requested a review from hidmic September 1, 2025 23:20

		template <class T, class U, class V>
		constexpr auto operator()(T&& t, U&& u, V& v) const {

		@@ -0,0 +1 @@
		test_2025-06-01_21-06-40 No newline at end of file

		/// 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,


		std::vector<double> lv_total_var, mult_total_var;

		for (int trial = 0; trial < std::min(50, num_trials); ++trial) { // Subset for performance

	const auto sample_size = std::max(static_cast<std::size_t>(1), container_size / 10); // Sample 10% of particles
	const auto sample_size = std::max(1UL, container_size / 10UL); // Sample 10% of particles

Conversation

larodriguez22 commented Apr 30, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Proposed changes

Type of change

Checklist

Additional comments

How to run it?

Uh oh!

hidmic left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

hidmic left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

glpuga left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

larodriguez22 commented Jun 3, 2025

Uh oh!

hidmic left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

hidmic left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

larodriguez22 commented Apr 30, 2025 •

edited

Loading