chen_method

//chen.h
#ifndef CHEN_H_
#define CHEN_H_
#include <vector>
#include <boost/random.hpp>
/// \file chen.h
/// \mainpage Chen discrete finite distribution generator algorithm
namespace chen
{
class Chen {
public:
  explicit Chen(const std::vector<double>& probs);
  /**
   * Generates new pseudo-random value
   * \return is 0..n-1
   */
  size_t next();
  /**
   * Return the already created value
   * \sa next()
   */
  size_t lookup()const;
  void setProbs(const std::vector<double>& p);
  void setSeed(int seed) const;
private:
  size_t _cur_value;
  std::vector<double> _p;
  std::vector<double> _s;
  size_t _param;
  std::vector<size_t> _r;
  boost::uniform_01<> _alpha;
  mutable boost::minstd_rand _gen;
};
}//namespace chen
#endif
//chen.cpp
#include "chen.h"
#include <cmath>
#include <numeric>

namespace chen {
Chen::Chen(const std::vector<double>&probs): _param(probs.size()+1)
{
  setProbs(probs);
  next();
}

size_t Chen::next()
{
  double alpha = _alpha(_gen);
  size_t j=(size_t)std::floor(_param*alpha);
  size_t i=_r[j];
  while(alpha>=_s[i])
    ++i;
  _cur_value=i;
  return _cur_value;
}

void Chen::setProbs(const std::vector<double>&probs)
{
  _p=probs;
  _s=_p;
  std::partial_sum(_s.begin(),_s.end(),_s.begin());//s_i=sum(p_1+...p_i)
  size_t i = 0;
  double t=0;
  _param = probs.size()+1;
  _r.resize(_param);
  double revparam=1./_param;
  for(size_t j=0; j<_param; ++j) {
    while(_s[i]<=t)
      ++i;
    _r[j] = i;
    t += revparam;
  }
}

void Chen::setSeed(int seed)const
{
  _gen.seed(seed);
}

size_t Chen::lookup() const
{
  return _cur_value;
}
} //namespace chen