00001 #ifndef SCITBX_RANDOM_H
00002 #define SCITBX_RANDOM_H
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include <scitbx/math/r3_rotation.h>
00015 #include <scitbx/array_family/shared.h>
00016 #include <boost/cstdint.hpp>
00017 #include <stdexcept>
00018
00019 namespace scitbx {
00020 namespace boost_random {
00021
00022
00023 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
00024 int s, UIntType b, int t, UIntType c, int l, UIntType val>
00025 class mersenne_twister
00026 {
00027 public:
00028 typedef UIntType result_type;
00029 static const int word_size = w;
00030 static const int state_size = n;
00031 static const int shift_size = m;
00032 static const int mask_bits = r;
00033 static const UIntType parameter_a = a;
00034 static const int output_u = u;
00035 static const int output_s = s;
00036 static const UIntType output_b = b;
00037 static const int output_t = t;
00038 static const UIntType output_c = c;
00039 static const int output_l = l;
00040
00041 static const bool has_fixed_range = false;
00042
00043 mersenne_twister() { seed(); }
00044
00045 explicit mersenne_twister(const UIntType& value)
00046 { seed(value); }
00047
00048 void seed() { seed(UIntType(5489)); }
00049
00050 void seed(const UIntType& value)
00051 {
00052
00053
00054
00055
00056 const UIntType mask = ~0u;
00057 x[0] = value & mask;
00058 for (i = 1; i < n; i++) {
00059
00060 x[i] = (1812433253UL * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask;
00061 }
00062 }
00063
00064 result_type min_value() const { return 0; }
00065 result_type max_value() const
00066 {
00067
00068 result_type res = 0;
00069 for(int j = 0; j < w; ++j)
00070 res |= (1u << j);
00071 return res;
00072 }
00073
00074 result_type operator()();
00075
00076 af::shared<std::size_t>
00077 getstate() const
00078 {
00079 af::shared<std::size_t> result;
00080 result.reserve(n);
00081 for(unsigned j=0;j<n;j++) {
00082 result.push_back(compute(j));
00083 }
00084 return result;
00085 }
00086
00087 void
00088 setstate(af::const_ref<std::size_t> const& state)
00089 {
00090 if (state.size() != n) {
00091 throw std::runtime_error(
00092 "mersenne_twister::setstate: improper state.size()");
00093 }
00094 for(unsigned j=0;j<n;j++) {
00095 x[j] = state[j];
00096 }
00097 i = n;
00098 }
00099
00100 private:
00101
00102 UIntType compute(unsigned int index) const
00103 {
00104
00105 return x[ (i + n + index) % (2*n) ];
00106 }
00107 void twist(int block);
00108
00109
00110
00111
00112
00113
00114
00115 UIntType x[2*n];
00116 int i;
00117 };
00118
00119 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
00120 int s, UIntType b, int t, UIntType c, int l, UIntType val>
00121 void mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::twist(int block)
00122 {
00123 const UIntType upper_mask = (~0u) << r;
00124 const UIntType lower_mask = ~upper_mask;
00125
00126 if(block == 0) {
00127 for(int j = n; j < 2*n; j++) {
00128 UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask);
00129 x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
00130 }
00131 } else if (block == 1) {
00132
00133 {
00134 for(int j = 0; j < n-m; j++) {
00135 UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
00136 x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0);
00137 }
00138 }
00139
00140 for(int j = n-m; j < n-1; j++) {
00141 UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
00142 x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
00143 }
00144
00145 UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask);
00146 x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0);
00147 i = 0;
00148 }
00149 }
00150
00151 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
00152 int s, UIntType b, int t, UIntType c, int l, UIntType val>
00153 inline typename mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::result_type
00154 mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::operator()()
00155 {
00156 if(i == n)
00157 twist(0);
00158 else if(i >= 2*n)
00159 twist(1);
00160
00161 UIntType z = x[i];
00162 ++i;
00163 z ^= (z >> u);
00164 z ^= ((z << s) & b);
00165 z ^= ((z << t) & c);
00166 z ^= (z >> l);
00167 return z;
00168 }
00169
00170
00171 typedef mersenne_twister<boost::uint32_t,32,624,397,31,0x9908b0df,11,
00172 7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937;
00173
00174 }
00175 }
00176
00177 namespace scitbx {
00179
00181 namespace random {
00182
00184
00186 class mersenne_twister
00187 {
00188 public:
00190 mersenne_twister(unsigned seed=0)
00191 :
00192 generator_(seed+1)
00193 {
00194
00195 }
00196
00198 void
00199 seed(unsigned value=0) { generator_.seed(value+1); }
00200
00202 std::size_t
00203 random_size_t_min()
00204 {
00205 return static_cast<std::size_t>(generator_.min_value());
00206 }
00207
00209 std::size_t
00210 random_size_t_max()
00211 {
00212 return static_cast<std::size_t>(generator_.max_value());
00213 }
00214
00218 std::size_t
00219 random_size_t()
00220 {
00221 return static_cast<std::size_t>(generator_());
00222 }
00223
00227 af::shared<std::size_t>
00228 random_size_t(std::size_t size)
00229 {
00230 af::shared<std::size_t> result(
00231 size, af::init_functor_null<std::size_t>());
00232 std::size_t* r = result.begin();
00233 for(std::size_t i=0;i<size;i++) *r++ = random_size_t();
00234 return result;
00235 }
00236
00240 af::shared<std::size_t>
00241 random_size_t(std::size_t size, std::size_t modulus)
00242 {
00243 af::shared<std::size_t> result(
00244 size, af::init_functor_null<std::size_t>());
00245 std::size_t* r = result.begin();
00246 for(std::size_t i=0;i<size;i++) *r++ = random_size_t() % modulus;
00247 return result;
00248 }
00249
00253 double
00254 random_double()
00255 {
00256
00257
00258
00259
00260
00261
00262
00263
00264 std::size_t a = random_size_t() >> 5;
00265 std::size_t b = random_size_t() >> 6;
00266 static const double c = 1.0/9007199254740992.0;
00267 return (a*67108864.0+b)*c;
00268 }
00269
00273 af::shared<double>
00274 random_double(std::size_t size)
00275 {
00276 af::shared<double> result(size, af::init_functor_null<double>());
00277 double* r = result.begin();
00278 for(std::size_t i=0;i<size;i++) *r++ = random_double();
00279 return result;
00280 }
00281
00285 af::shared<double>
00286 random_double(std::size_t size, double factor)
00287 {
00288 af::shared<double> result(size, af::init_functor_null<double>());
00289 double* r = result.begin();
00290 for(std::size_t i=0;i<size;i++) *r++ = random_double() * factor;
00291 return result;
00292 }
00293
00295 af::shared<bool>
00296 random_bool(std::size_t size, double threshold)
00297 {
00298 af::shared<bool> result(size, af::init_functor_null<bool>());
00299 bool *r = result.begin();
00300 for(std::size_t i=0;i<size;i++) *r++ = random_double() < threshold;
00301 return result;
00302 }
00303
00305 af::shared<std::size_t>
00306 random_permutation(std::size_t size)
00307 {
00308 af::shared<std::size_t> result(
00309 size, af::init_functor_null<std::size_t>());
00310 std::size_t* r = result.begin();
00311 for(std::size_t i=0;i<size;i++) *r++ = i;
00312 r = result.begin();
00313 for(std::size_t i=0;i<size;i++) {
00314 std::size_t j = static_cast<std::size_t>(generator_()) % size;
00315 std::swap(r[i], r[j]);
00316 }
00317 return result;
00318 }
00319
00321
00340 scitbx::vec3<double>
00341 random_double_point_on_sphere()
00342 {
00343 vec3<double> result;
00344 double z = 2 * random_double() - 1;
00345 double t = constants::two_pi * random_double();
00346 double r = std::sqrt(1-z*z);
00347 result[0] = r * std::cos(t);
00348 result[1] = r * std::sin(t);
00349 result[2] = z;
00350 return result;
00351 }
00352
00353 af::tiny<double, 4>
00354 random_double_unit_quaternion()
00355 {
00356
00357
00358
00359
00360
00361
00362
00363
00364 scitbx::vec3<double> axis = random_double_point_on_sphere();
00365 double angle = constants::two_pi * random_double();
00366 return math::r3_rotation::axis_and_angle_as_unit_quaternion(
00367 axis, angle, false);
00368 }
00369
00370 scitbx::mat3<double>
00371 random_double_r3_rotation_matrix()
00372 {
00373
00374 scitbx::vec3<double> axis = random_double_point_on_sphere();
00375 double angle = constants::two_pi * random_double();
00376 return math::r3_rotation::axis_and_angle_as_matrix(
00377 axis, angle, false);
00378 }
00379
00380 af::shared<std::size_t>
00381 getstate() const { return generator_.getstate(); }
00382
00383 void
00384 setstate(af::const_ref<std::size_t> const& state)
00385 {
00386 generator_.setstate(state);
00387 }
00388
00389 private:
00390 boost_random::mt19937 generator_;
00391 };
00392
00393 }}
00394
00395 #endif // SCITBX_RANDOM_H