<commit_before>#ifndef VSMC_CORE_RESAMPLE_HPP
#define VSMC_CORE_RESAMPLE_HPP

#include <vsmc/internal/common.hpp>

namespace vsmc {

namespace internal {

template <typename SizeType, typename RngSetType>
inline void weight2replication (SizeType N, SizeType S, RngSetType &rng_set,
        const double *weight, SizeType *replication)
{
        double sum_w = 0;
        double acc_w = 0;
        SizeType acc_s = 0;

        for (SizeType i = 0; i != N; ++i) {
            replication[i] = 0;
            sum_w += weight[i];
        }

        for (SizeType i = 0; i != N; ++i) {
            if (acc_s < S && weight[i] > 0) {
                typedef typename cxx11::make_signed<SizeType>::type s_t;
                s_t s = S - acc_s;
                double p = weight[i] / (sum_w - acc_w);
                if (p < 0) {
                    p = 0;
                    assert(p > -1e-6);
                }
                if (p > 1) {
                    p = 1;
                    assert(p - 1 < 1e-6);
                }
                cxx11::binomial_distribution<s_t> binom(s, p);
                replication[i] = binom(rng_set.rng(i));
            }
            acc_w += weight[i];
            acc_s += replication[i];
        }
}

} // namespace vsmc::internal

/// \brief Resample scheme
/// \ingroup Core
enum ResampleScheme {
    MULTINOMIAL,         ///< Multinomial resampling
    RESIDUAL,            ///< Reisudal resampling
    STRATIFIED,          ///< Startified resampling
    SYSTEMATIC,          ///< Systematic resampling
    RESIDUAL_STRATIFIED, ///< Stratified resampling on the residuals
    RESIDUAL_SYSTEMATIC  ///< Systematic resampling on the residuals
}; // enum ResamleScheme

/// \brief Int-to-Type struct template for resampling scheme
/// \ingroup Core
template <typename EnumType, EnumType S> struct ResampleType {};

/// \brief Resample class template
/// \ingroup Core
template <typename ResType, typename SizeType, typename RngSetType>
class Resample {};

/// \brief Multinomial resampling
/// \ingroup Core
template <typename SizeType, typename RngSetType>
class Resample<ResampleType<ResampleScheme, MULTINOMIAL>,
      SizeType, RngSetType>
{
    public :

    void operator() (SizeType N, RngSetType &rng_set,
            const double *weight, SizeType *replication)
    {
        internal::weight2replication(N, N, rng_set, weight, replication);
    }
}; // Mulitnomial resampling

/// \brief Residual resampling
/// \ingroup Core
template <typename SizeType, typename RngSetType>
class Resample<ResampleType<ResampleScheme, RESIDUAL>,
      SizeType, RngSetType>
{
    public :

    void operator() (SizeType N, RngSetType &rng_set,
            const double *weight, SizeType *replication)
    {
        using std::modf;

        residual_.resize(N);
        integral_.resize(N);
        for (SizeType i = 0; i != N; ++i)
            residual_[i] = modf(N * weight[i], &integral_[i]);
        SizeType S = static_cast<SizeType>(residual_.sum());
        internal::weight2replication(N, S, rng_set, weight, replication);<commit_msg>fix residual resampling<commit_after>#ifndef VSMC_CORE_RESAMPLE_HPP
#define VSMC_CORE_RESAMPLE_HPP

#include <vsmc/internal/common.hpp>

namespace vsmc {

namespace internal {

template <typename SizeType, typename RngSetType>
inline void weight2replication (SizeType N, SizeType S, RngSetType &rng_set,
        const double *weight, SizeType *replication)
{
        double sum_w = 0;
        double acc_w = 0;
        SizeType acc_s = 0;

        for (SizeType i = 0; i != N; ++i) {
            replication[i] = 0;
            sum_w += weight[i];
        }

        for (SizeType i = 0; i != N; ++i) {
            if (acc_s < S && weight[i] > 0) {
                typedef typename cxx11::make_signed<SizeType>::type s_t;
                s_t s = S - acc_s;
                double p = weight[i] / (sum_w - acc_w);
                if (p < 0) {
                    p = 0;
                    assert(p > -1e-6);
                }
                if (p > 1) {
                    p = 1;
                    assert(p - 1 < 1e-6);
                }
                cxx11::binomial_distribution<s_t> binom(s, p);
                replication[i] = binom(rng_set.rng(i));
            }
            acc_w += weight[i];
            acc_s += replication[i];
        }
}

} // namespace vsmc::internal

/// \brief Resample scheme
/// \ingroup Core
enum ResampleScheme {
    MULTINOMIAL,         ///< Multinomial resampling
    RESIDUAL,            ///< Reisudal resampling
    STRATIFIED,          ///< Startified resampling
    SYSTEMATIC,          ///< Systematic resampling
    RESIDUAL_STRATIFIED, ///< Stratified resampling on the residuals
    RESIDUAL_SYSTEMATIC  ///< Systematic resampling on the residuals
}; // enum ResamleScheme

/// \brief Int-to-Type struct template for resampling scheme
/// \ingroup Core
template <typename EnumType, EnumType S> struct ResampleType {};

/// \brief Resample class template
/// \ingroup Core
template <typename ResType, typename SizeType, typename RngSetType>
class Resample {};

/// \brief Multinomial resampling
/// \ingroup Core
template <typename SizeType, typename RngSetType>
class Resample<ResampleType<ResampleScheme, MULTINOMIAL>,
      SizeType, RngSetType>
{
    public :

    void operator() (SizeType N, RngSetType &rng_set,
            const double *weight, SizeType *replication)
    {
        internal::weight2replication(N, N, rng_set, weight, replication);
    }
}; // Mulitnomial resampling

/// \brief Residual resampling
/// \ingroup Core
template <typename SizeType, typename RngSetType>
class Resample<ResampleType<ResampleScheme, RESIDUAL>,
      SizeType, RngSetType>
{
    public :

    void operator() (SizeType N, RngSetType &rng_set,
            const double *weight, SizeType *replication)
    {
        using std::modf;

        residual_.resize(N);
        integral_.resize(N);
        for (SizeType i = 0; i != N; ++i)
            residual_[i] = modf(N * weight[i], &integral_[i]);
        double dsize = residual_.sum();
        SizeType size = static_cast<SizeType>(dsize);
        residual_ /= dsize;
        internal::weight2replication(N, size, rng_set, residual_.data(),
                replication);
{'fraction_non_alphanumeric': 0.11342518423582185, 'fraction_numerical': 0.00768984299903877, 'mean_word_length': 2.1340361445783134, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4528989', 'n_tokens_mistral': 2064, 'n_tokens_neox': 1943, 'n_words': 502}