1 ////////////////////////////////////////////////////////////////////////////////////// 2 // This file is distributed under the University of Illinois/NCSA Open Source License. 3 // See LICENSE file in top directory for details. 4 // 5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. 6 // 7 // File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 8 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign 9 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory 10 // 11 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 12 ////////////////////////////////////////////////////////////////////////////////////// 13 14 15 #ifndef OHMMS_BOOSTRANDOM_H 16 #define OHMMS_BOOSTRANDOM_H 17 18 #ifdef HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 #include <ctime> 22 #include <sstream> 23 #include <limits> 24 #include <boost/config.hpp> 25 #ifdef BOOST_CLANG 26 #pragma clang diagnostic push 27 #pragma clang diagnostic ignored "-W#pragma-messages" 28 #endif 29 #include <boost/random.hpp> 30 #ifdef BOOST_CLANG 31 #pragma clang diagnostic pop 32 #endif 33 34 /** random number generator using boost::random 35 * 36 * A wrapper of boost::random class to work with applicatoins. 37 */ 38 template<typename T, typename RNG = boost::mt19937> 39 class BoostRandom 40 { 41 public: 42 /// real result type 43 typedef T result_type; 44 /// randmon number generator [0,max) where max depends on the generator type 45 typedef RNG generator_type; 46 /// unsigned integer type 47 typedef typename generator_type::result_type uint_type; 48 /// real random generator [0,1) 49 typedef boost::variate_generator<generator_type, boost::uniform_real<T>> uniform_generator_type; 50 51 std::string ClassName; 52 std::string EngineName; 53 54 ///default constructor 55 explicit BoostRandom(uint_type iseed = 911, const std::string& aname = "mt19937") 56 : ClassName("boost"), 57 EngineName(aname), 58 myContext(0), 59 nContexts(1), 60 baseOffset(0), 61 uni(generator_type(iseed), boost::uniform_real<T>(0, 1)) 62 {} 63 64 ///copy constructor BoostRandom(const BoostRandom & rng)65 BoostRandom(const BoostRandom& rng) 66 : ClassName(rng.ClassName), 67 EngineName(rng.EngineName), 68 myContext(rng.myContext), 69 nContexts(rng.nContexts), 70 baseOffset(rng.baseOffset), 71 uni(rng.uni) 72 {} 73 74 ///copy operator (unnecessary but why not) 75 BoostRandom<T, RNG>& operator=(const BoostRandom& r) 76 { 77 ClassName = r.ClassName; 78 EngineName = r.EngineName; 79 myContext = r.myContext; 80 nContexts = r.nContexts; 81 baseOffset = r.baseOffset, uni = r.uni; 82 return *this; 83 } 84 ~BoostRandom()85 ~BoostRandom() {} 86 87 /** initialize the generator 88 * @param i thread index 89 * @param nstr number of threads 90 * @param iseed_in input seed 91 * 92 * Initialize generator with the seed. 93 */ 94 void init(int i, int nstr, int iseed_in, uint_type offset = 1) 95 { 96 uint_type baseSeed = iseed_in; 97 myContext = i; 98 nContexts = nstr; 99 if (iseed_in <= 0) 100 baseSeed = make_seed(i, nstr); 101 baseOffset = offset; 102 uni.engine().seed(baseSeed); 103 } 104 105 ///get baseOffset offset()106 inline int offset() const { return baseOffset; } 107 ///assign baseOffset offset()108 inline int& offset() { return baseOffset; } 109 110 ///assign seed seed(uint_type aseed)111 inline void seed(uint_type aseed) { uni.engine().seed(aseed); } 112 engine()113 uniform_generator_type& engine() { return uni; } 114 /////reset the seed 115 //inline void reset() 116 //{ 117 // uni.engine().seed(make_seed(myContext,nContexts,baseOffset)); 118 //} 119 120 /** return a random number [0,1) 121 */ rand()122 inline result_type rand() { return uni(); } 123 124 /** return a random number [0,1) 125 */ operator()126 inline result_type operator()() { return uni(); } 127 128 /** return a random integer 129 */ irand()130 inline uint_type irand() { return uni.engine()() % std::numeric_limits<uint_type>::max(); } 131 132 /** generate a series of random numbers */ 133 template<typename T1> generate_uniform(T1 * restrict d,int n)134 inline void generate_uniform(T1* restrict d, int n) 135 { 136 for (int i = 0; i < n; ++i) 137 d[i] = uni(); 138 } 139 generate_normal(T * restrict d,int n)140 inline void generate_normal(T* restrict d, int n) { BoxMuller2::generate(*this, d, n); } 141 142 //inline void bivariate(resul_type& g1, resul_type &g2) { 143 // resul_type v1, v2, r; 144 // do { 145 // v1 = 2.0e0*uni() - 1.0e0; 146 // v2 = 2.0e0*uni() - 1.0e0; 147 // r = v1*v1+v2*v2; 148 // } while(r > 1.0e0); 149 // resul_type fac = sqrt(-2.0e0*log(r)/r); 150 // g1 = v1*fac; 151 // g2 = v2*fac; 152 //} 153 state_size()154 inline size_t state_size() const { return uni.engine().state_size; } 155 read(std::istream & rin)156 inline void read(std::istream& rin) { rin >> uni.engine(); } 157 write(std::ostream & rout)158 inline void write(std::ostream& rout) const { rout << uni.engine(); } 159 save(std::vector<uint_type> & curstate)160 inline void save(std::vector<uint_type>& curstate) const 161 { 162 curstate.clear(); 163 std::stringstream otemp; 164 otemp << uni.engine(); 165 copy(std::istream_iterator<uint_type>(otemp), std::istream_iterator<uint_type>(), back_inserter(curstate)); 166 } 167 load(const std::vector<uint_type> & newstate)168 inline void load(const std::vector<uint_type>& newstate) 169 { 170 std::stringstream otemp; 171 copy(newstate.begin(), newstate.end(), std::ostream_iterator<uint_type>(otemp, " ")); 172 otemp >> uni.engine(); 173 } 174 175 private: 176 ///context number 177 int myContext; 178 ///number of contexts 179 int nContexts; 180 ///offset of the random seed 181 int baseOffset; 182 ///random number generator [0,1) 183 uniform_generator_type uni; 184 }; 185 #endif 186