sampling.h

   1 #ifndef SAMPLING
   2 #define SAMPLING
   3
   4 #include<ctime>
   5 #include<cstdio>
   6 #include<cassert>
   7 #include<vector>
   8 #include<set>
   9
  10 #include "boost/random.hpp"
  11
  12 typedef unsigned int seedType;
  13 typedef boost::mt19937 engine_type;
  14 typedef boost::gamma_distribution<> gamma_dist;
  15 typedef boost::uniform_01<engine_type> uniform01;
  16 typedef boost::variate_generator<engine_type&, gamma_dist> gamma_generator;
  17
  18 class engineFactory {
  19 public:
  20   static void init() { seedEngine = new engine_type(time(NULL)); }
  21   static void init(seedType seed) { seedEngine = new engine_type(seed); }
  22
  23   static void finish() { if (seedEngine != NULL) delete seedEngine; }
  24
  25         static engine_type *new_engine() {
  26                 seedType seed;
  27                 static std::set<seedType> seedSet;                      // empty set of seeds
  28                 std::set<seedType>::iterator iter;
  29
  30                 do {
  31                         seed = (*seedEngine)();
  32                         iter = seedSet.find(seed);
  33                 } while (iter != seedSet.end());
  34                 seedSet.insert(seed);
  35
  36                 return new engine_type(seed);
  37         }
  38
  39  private:
  40         static engine_type *seedEngine;
  41 };
  42
  43 engine_type* engineFactory::seedEngine = NULL;
  44
  45 // arr should be cumulative!
  46 // interval : [,)
  47 // random number should be in [0, arr[len - 1])
  48 // If by chance arr[len - 1] == 0.0, one possibility is to sample uniformly from 0...len-1
  49 int sample(uniform01& rg, std::vector<double>& arr, int len) {
  50   int l, r, mid;
  51   double prb = rg() * arr[len - 1];
  52
  53   l = 0; r = len - 1;
  54   while (l <= r) {
  55     mid = (l + r) / 2;
  56     if (arr[mid] <= prb) l = mid + 1;
  57     else r = mid - 1;
  58   }
  59
  60   if (l >= len) { printf("%d %lf %lf\n", len, arr[len - 1], prb); }
  61   assert(l < len);
  62
  63   return l;
  64 }
  65
  66 #endif