void ebpps_sample::merge()

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_));
    }
  }
}