RngThread.H

Go to the documentation of this file.

 
#ifndef RNGTHREAD_H_
#define RNGTHREAD_H_
 
#include <random>
#include <sstream>
#include <vector>
#include "ThreadVec.H"
#include "Types.H"
#include "../pcg-cpp/pcg_random.hpp"
#include "../IO/ParmParser.H"
#ifdef _OPENMP
#   include <omp.h>
#endif
 
namespace criptic {
 
  class RngThread {
 
  public:
 
    RngThread(const ParmParser &pp);
 
    RngThread(std::string& buf);
 
    ~RngThread();
 
    Real uniform(const Real a = 0.0, const Real b = 1.0) const {
      std::uniform_real_distribution<Real> dist(a,b);
      pcg64& rth = *r();
      return dist(rth);
    }
 
    IdxType uniformInt(const IdxType a = 0, const IdxType b = maxIdx) const {
      std::uniform_int_distribution<IdxType> dist(a,b);
      pcg64& rth = *r();
      return dist(rth);
    }
    
    Real normal(const Real mean = 0.0, const Real stddev = 1.0) const {
      std::normal_distribution<Real> dist(mean, stddev);
      pcg64& rth = *r();
      return dist(rth);
    }
 
    IdxType poisson(const Real ex) const {
      std::poisson_distribution<IdxType> dist(ex);
      pcg64& rth = *r();
      return dist(rth);
    }
 
    template <typename T>
    T& choose(std::vector<T>& elem,
          const std::vector<Real>& w = std::vector<Real>()) const {
      if (elem.size() == 1) return elem[0];
      else if (w.size() == 0) return elem[uniformInt(0, elem.size()-1)];
      else {
    std::vector<Real> wCum(w.size());
    wCum[0] = w[0];
    for (IdxType i = 1; i < w.size(); i++) wCum[i] = wCum[i-1] + w[i];
    Real x = uniform(0, wCum.back());
    if (x < wCum[0]) return elem[0];
    else {
      IdxType l = 0, r = wCum.size()-1;
      while (r - l > 1) {
        IdxType c = (l + r) / 2;
        if (x < wCum[c]) r = c;
        else l = c;
      }
      return elem[r];
    }
      }
    }
    
    template <typename T>
    const T& choose(const std::vector<T>& elem,
          const std::vector<Real>& w = std::vector<Real>()) const {
      if (w.size() == 0) return elem[uniformInt(0, elem.size()-1)];
      else {
    std::vector<Real> wCum(w.size());
    wCum[0] = w[0];
    for (IdxType i = 0; i < w.size(); i++) wCum[i] = wCum[i-1] + w[i];
    Real x = uniform(0, wCum.back());
    if (x < wCum[0]) return elem[0];
    else {
      IdxType l = 0, r = wCum.size()-1;
      while (r - l > 1) {
        IdxType c = (l + r) / 2;
        if (x < wCum[c]) r = c;
        else l = c;
      }
      return elem[r];
    }
      }
    }
 
    const std::string getState() const;
 
    void setState(const std::string& buf);
 
    static IdxType stateSize(); 
 
  private:
 
    // Internal data
    ThreadVec<pcg64 *> r;  
  };
 
}
 
#endif
// _RNGTHREAD_H_