crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-21 06:25:30 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

rejection sampling for lin-log interpolation

Browse source
This commit is contained in:
myerspat 2023-05-18 15:15:25 -04:00
parent b98c4590d7
commit 7ba6a45001
1 changed files with 22 additions and 1 deletions
Download patch file Download diff file Expand all files Collapse all files

23
src/distribution.cpp
View file

@ -304,7 +304,7 @@ void Tabular::init(
c_[i] = c_[i - 1] + p_[i - 1] * (x_[i] - x_[i - 1]) /
std::log(p_[i] / p_[i - 1]) *
(p_[i] / p_[i - 1] - 1);
} else {
} else if (interp_ == Interpolation::log_log) {
double m = std::log(p_[i] / p_[i - 1] - x_[i] / x_[i - 1]);
c_[i] =
c_[i - 1] + p_[i - 1] / (std::pow(x_[i - 1], m) * (m + 1)) *
@ -359,6 +359,27 @@ double Tabular::sample(uint64_t* seed) const
(std::sqrt(std::max(0.0, p_i * p_i + 2 * m * (c - c_i))) - p_i) /
m;
}
} else if (interp_ == Interpolation::lin_log) {
// Linear-log interpolation
// Inverse transform sampling with linear-log interpolation is not possible
// so this uses rejection sampling
double x_i1 = x_[i + 1];
double p_i1 = p_[i + 1];
if (p_i1 == p_i) {
return x_i + (c - c_i) / p_i;
}
while (true) {
double x = x_i + prn(seed) * (x_i1 - x_i);
double p_unif = p_i + prn(seed) * (p_i1 - p_i);
// Sample linear-log PDF
double p = p_i + std::log(x / x_i - x_i1 / x_i) * (p_i1 - p_i);
if (p_unif < p) {
return x;
}
}
} else if (interp_ == Interpolation::log_lin) {
// Log-linear interpolation
double x_i1 = x_[i + 1];