sampling.h

#ifndef SAMPLING
#define SAMPLING

#include <ctime>
#include <cstdio>
#include <cassert>
#include <vector>
#include <set>
#include <random>

#include "boost_compat.h"

typedef unsigned int seedType;
typedef std::mt19937 engine_type;

// Boost-compatible distributions (match Boost.Random for reproducible results)
struct uniform_01_dist {
	typedef double result_type;
	template<typename Engine>
	result_type operator()(Engine& eng) const { return boost_uniform_01(eng); }
};

struct gamma_dist {
	typedef double result_type;
	double alpha_, beta_;
	explicit gamma_dist(double alpha = 1.0, double beta = 1.0) : alpha_(alpha), beta_(beta) {}
	template<typename Engine>
	result_type operator()(Engine& eng) const { return boost_gamma(eng, alpha_, beta_); }
};

// Wrapper to provide variate_generator-like interface (callable with operator())
template<typename Engine, typename Dist>
class variate_generator {
	Engine& engine;
	Dist dist;
public:
	variate_generator(Engine& e, const Dist& d) : engine(e), dist(d) {}
	typename Dist::result_type operator()() { return dist(engine); }
};

typedef variate_generator<engine_type, uniform_01_dist> uniform_01_generator;
typedef variate_generator<engine_type, gamma_dist> gamma_generator;

class engineFactory {
public:
  static void init() { seedEngine = new engine_type(static_cast<engine_type::result_type>(time(NULL))); }
  static void init(seedType seed) { seedEngine = new engine_type(seed); }

  static void finish() { if (seedEngine != NULL) delete seedEngine; }

	static engine_type *new_engine() {
		seedType seed;
		static std::set<seedType> seedSet;			// empty set of seeds
		std::set<seedType>::iterator iter;

		do {
			seed = (*seedEngine)();
			iter = seedSet.find(seed);
		} while (iter != seedSet.end());
		seedSet.insert(seed);

		return new engine_type(seed);
	}

 private:
	static engine_type *seedEngine;
};

engine_type* engineFactory::seedEngine = NULL;

// arr should be cumulative!
// interval : [,)
// random number should be in [0, arr[len - 1])
// If by chance arr[len - 1] == 0.0, one possibility is to sample uniformly from 0...len-1
int sample(uniform_01_generator& rg, std::vector<double>& arr, int len) {
  int l, r, mid;
  double prb = rg() * arr[len - 1];

  l = 0; r = len - 1;
  while (l <= r) {
    mid = (l + r) / 2;
    if (arr[mid] <= prb) l = mid + 1;
    else r = mid - 1;
  }

  if (l >= len) { printf("%d %lf %lf\n", len, arr[len - 1], prb); }
  assert(l < len);

  return l;
}

#endif