in sampling/include/ebpps_sample_impl.hpp [123:169]
void ebpps_sample<T,A>::merge(FwdSample&& other) {
double c_int;
const double c_frac = std::modf(c_, &c_int);
double unused;
const double other_c_frac = std::modf(other.c_, &unused);
// update c_ here but do NOT recompute fractional part yet
c_ += other.c_;
for (uint32_t i = 0; i < other.data_.size(); ++i)
data_.emplace_back(conditional_forward<FwdSample>(other.data_[i]));
// This modifies the original algorithm slightly due to numeric
// precision issues. Specifically, the test if c_frac + other_c_frac == 1.0
// happens before tests for < 1.0 or > 1.0 and can also be triggered
// if c_ == floor(c_) (the updated value of c_, not the input).
//
// We can still run into issues where c_frac + other_c_frac == epsilon
// and the first case would have ideally triggered. As a result, we must
// check if the partial item exists before adding to the data_ vector.
if (c_frac == 0.0 && other_c_frac == 0.0) {
partial_item_.reset();
} else if (c_frac + other_c_frac == 1.0 || c_ == std::floor(c_)) {
if (next_double() <= c_frac) {
if (partial_item_)
data_.emplace_back(std::move(*partial_item_));
} else {
if (other.partial_item_)
data_.emplace_back(conditional_forward<FwdSample>(*other.partial_item_));
}
partial_item_.reset();
} else if (c_frac + other_c_frac < 1.0) {
if (next_double() > c_frac / (c_frac + other_c_frac)) {
set_partial(conditional_forward<FwdSample>(*other.partial_item_));
}
} else { // c_frac + other_c_frac > 1
if (next_double() <= (1 - c_frac) / ((1 - c_frac) + (1 - other_c_frac))) {
data_.emplace_back(conditional_forward<FwdSample>(*other.partial_item_));
} else {
data_.emplace_back(std::move(*partial_item_));
partial_item_.reset();
set_partial(conditional_forward<FwdSample>(*other.partial_item_));
}
}
}