crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Use binary search for tabular distribution sampling

Browse source
This commit is contained in:
Paul Romano 2026-07-14 18:57:16 -05:00
parent 23572ccc92
commit 7c1af32c5e
2 changed files with 28 additions and 8 deletions
Download patch file Download diff file Expand all files Collapse all files

11
src/distribution.cpp
View file

@ -552,14 +552,9 @@ double Tabular::sample_unbiased(uint64_t* seed) const
double c = prn(seed);
// Find first CDF bin which is above the sampled value
double c_i = c_[0];
int i;
std::size_t n = c_.size();
for (i = 0; i < n - 1; ++i) {
if (c <= c_[i + 1])
break;
c_i = c_[i + 1];
}
auto c_iter = std::lower_bound(c_.begin() + 1, c_.end(), c);
int i = std::distance(c_.begin(), c_iter) - 1;
double c_i = c_[i];
// Determine bounding PDF values
double x_i = x_[i];

25
tests/cpp_unit_tests/test_distribution.cpp
View file

@ -82,6 +82,31 @@ TEST_CASE("Test alias sampling method for pugixml constructor")
}
}
TEST_CASE("Test sampling a large linear-linear tabular distribution")
{
constexpr int n_points = 10001;
constexpr int n_samples = 200000;
openmc::vector<double> x(n_points);
openmc::vector<double> p(n_points);
for (int i = 0; i < n_points; ++i) {
x[i] = static_cast<double>(i) / (n_points - 1);
p[i] = 2.0 * x[i];
}
openmc::Tabular dist(
x.data(), p.data(), n_points, openmc::Interpolation::lin_lin);
uint64_t seed = openmc::init_seed(0, 0);
double mean = 0.0;
for (int i = 0; i < n_samples; ++i) {
mean += dist.sample(&seed).first;
}
mean /= n_samples;
// The normalized PDF is 2x on [0, 1], which has a mean of 2/3.
REQUIRE_THAT(mean, Catch::Matchers::WithinAbs(2.0 / 3.0, 0.003));
}
TEST_CASE("Test construction of SpatialBox with parameters")
{
openmc::Position ll {-1, -2, -3};