1 // -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
2 //
3 // Copyright (C) 2013-2014 Conrad Sanderson
4 // Copyright (C) 2013-2014 NICTA (www.nicta.com.au)
5 //
6 // This Source Code Form is subject to the terms of the Mozilla Public
7 // License, v. 2.0. If a copy of the MPL was not distributed with this
8 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 //
10 // This file is based on Conrad's default generators and as such licensed under both
11 // the MPL 2.0 for his as well as the GNU GPL 2.0 or later for my modification to it.
12
13 // Copyright (C) 2014 Dirk Eddelbuettel
14 //
15 // This file is part of RcppArmadillo.
16 //
17 // RcppArmadillo is free software: you can redistribute it and/or modify it
18 // under the terms of the GNU General Public License as published by
19 // the Free Software Foundation, either version 2 of the License, or
20 // (at your option) any later version.
21 //
22 // RcppArmadillo is distributed in the hope that it will be useful, but
23 // WITHOUT ANY WARRANTY; without even the implied warranty of
24 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 // GNU General Public License for more details.
26 //
27 // You should have received a copy of the GNU General Public License
28 // along with RcppArmadillo. If not, see <http://www.gnu.org/licenses/>.
29
30
31 // NB This files use R's uniform generator and can be compiled only when the R
32 // headers are available as is the case for RcppArmadillo.
33 //
34 // Also note that you MUST set / reset the R RNG state. When using RcppArmadillo
35 // via Rcpp Atttributes or the inline package, the RNGScope object is added which
36 // ensure this automatically. Should you build by hand, and omit both RNGScope as
37 // as manual calls to GetRNGState() and PutRNGState() you may get unstable results.
38 //
39 // See http://cran.r-project.org/doc/manuals/r-devel/R-exts.html#Random-numbers
40
41 class arma_rng_alt {
42 public:
43
44 typedef unsigned int seed_type;
45
46 inline static void set_seed(const seed_type val);
47
48 arma_inline static int randi_val();
49 arma_inline static double randu_val();
50 inline static double randn_val();
51
52 template<typename eT>
53 inline static void randn_dual_val(eT& out1, eT& out2);
54
55 template<typename eT>
56 inline static void randi_fill(eT* mem, const uword N, const int a, const int b);
57
58 inline static int randi_max_val();
59 };
60
set_seed(const arma_rng_alt::seed_type val)61 inline void arma_rng_alt::set_seed(const arma_rng_alt::seed_type val) {
62 // null-op, cannot set seed in R from C level code
63 // see http://cran.r-project.org/doc/manuals/r-devel/R-exts.html#Random-numbers
64 //
65 // std::srand(val);
66 (void) val; // to suppress a -Wunused warning
67 //
68 static int havewarned = 0;
69 if (havewarned++ == 0) {
70 ::Rf_warning("When called from R, the RNG seed has to be set at the R level via set.seed()");
71 }
72 }
73
randi_val()74 arma_inline int arma_rng_alt::randi_val() {
75 return ::Rf_runif(0, RAND_MAX); //std::rand();
76 }
77
randu_val()78 arma_inline double arma_rng_alt::randu_val() {
79 return double(::Rf_runif(0, 1));
80 //return double( double(std::rand()) * ( double(1) / double(RAND_MAX) ) );
81 }
82
randn_val()83 inline double arma_rng_alt::randn_val() {
84 // polar form of the Box-Muller transformation:
85 // http://en.wikipedia.org/wiki/Box-Muller_transformation
86 // http://en.wikipedia.org/wiki/Marsaglia_polar_method
87
88 double tmp1;
89 double tmp2;
90 double w;
91
92 do {
93 tmp1 = double(2) * double(::Rf_runif(0, 1)) - double(1);
94 tmp2 = double(2) * double(::Rf_runif(0, 1)) - double(1);
95 //tmp1 = double(2) * double(std::rand()) * (double(1) / double(RAND_MAX)) - double(1);
96 //tmp2 = double(2) * double(std::rand()) * (double(1) / double(RAND_MAX)) - double(1);
97
98 w = tmp1*tmp1 + tmp2*tmp2;
99 } while ( w >= double(1) );
100
101 return double( tmp1 * std::sqrt( (double(-2) * std::log(w)) / w) );
102 }
103
104 template<typename eT>
randn_dual_val(eT & out1,eT & out2)105 inline void arma_rng_alt::randn_dual_val(eT& out1, eT& out2) {
106 // make sure we are internally using at least floats
107 typedef typename promote_type<eT,float>::result eTp;
108
109 eTp tmp1;
110 eTp tmp2;
111 eTp w;
112
113 do {
114 tmp1 = eTp(2) * eTp(::Rf_runif(0, RAND_MAX)) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
115 tmp2 = eTp(2) * eTp(::Rf_runif(0, RAND_MAX)) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
116 //tmp1 = eTp(2) * eTp(std::rand()) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
117 //tmp2 = eTp(2) * eTp(std::rand()) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
118
119 w = tmp1*tmp1 + tmp2*tmp2;
120 } while ( w >= eTp(1) );
121
122 const eTp k = std::sqrt( (eTp(-2) * std::log(w)) / w);
123
124 out1 = eT(tmp1*k);
125 out2 = eT(tmp2*k);
126 }
127
128
129
130 template<typename eT>
randi_fill(eT * mem,const uword N,const int a,const int b)131 inline void arma_rng_alt::randi_fill(eT* mem, const uword N, const int a, const int b) {
132 if( (a == 0) && (b == RAND_MAX) ) {
133 for(uword i=0; i<N; ++i) {
134 mem[i] = ::Rf_runif(0, RAND_MAX);
135 //mem[i] = std::rand();
136 }
137 } else {
138 const uword length = b - a + 1;
139
140 const double scale = double(length) / double(RAND_MAX);
141
142 for(uword i=0; i<N; ++i) {
143 mem[i] = (std::min)( b, (int( double(::Rf_runif(0,RAND_MAX)) * scale ) + a) );
144 //mem[i] = (std::min)( b, (int( double(std::rand()) * scale ) + a) );
145 }
146 }
147 }
148
randi_max_val()149 inline int arma_rng_alt::randi_max_val() {
150 return RAND_MAX;
151 }
152