1
- // -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
2
- //
1
+
3
2
// Copyright (C) 2013-2014 Conrad Sanderson
4
3
// Copyright (C) 2013-2014 NICTA (www.nicta.com.au)
5
- //
4
+ //
6
5
// This Source Code Form is subject to the terms of the Mozilla Public
7
6
// License, v. 2.0. If a copy of the MPL was not distributed with this
8
7
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
8
//
10
- // This file is based on Conrad's default generators and as such licensed under both
9
+ // This file is based on Conrad's default generators and as such licensed under both
11
10
// the MPL 2.0 for his as well as the GNU GPL 2.0 or later for my modification to it.
12
11
13
12
// Copyright (C) 2014 Dirk Eddelbuettel
28
27
// along with RcppArmadillo. If not, see <http://www.gnu.org/licenses/>.
29
28
30
29
31
- // NB This files use R's uniform generator and can be compiled only when the R
30
+ // NB This files use R's uniform generator and can be compiled only when the R
32
31
// headers are available as is the case for RcppArmadillo.
33
32
//
34
33
// Also note that you MUST set / reset the R RNG state. When using RcppArmadillo
40
39
41
40
class arma_rng_alt {
42
41
public:
43
-
42
+
44
43
typedef unsigned int seed_type;
45
-
44
+
46
45
inline static void set_seed (const seed_type val);
47
-
46
+
48
47
arma_inline static int randi_val ();
49
48
arma_inline static double randu_val ();
50
49
inline static double randn_val ();
51
-
50
+
52
51
template <typename eT>
53
52
inline static void randn_dual_val (eT& out1, eT& out2);
54
-
53
+
55
54
template <typename eT>
56
55
inline static void randi_fill (eT* mem, const uword N, const int a, const int b);
57
-
56
+
58
57
inline static int randi_max_val ();
59
58
};
60
59
@@ -72,7 +71,7 @@ inline void arma_rng_alt::set_seed(const arma_rng_alt::seed_type val) {
72
71
}
73
72
74
73
arma_inline int arma_rng_alt::randi_val () {
75
- return ::Rf_runif (0 , RAND_MAX); // std::rand();
74
+ return static_cast < int >( ::Rf_runif (0 , RAND_MAX) ); // std::rand();
76
75
}
77
76
78
77
arma_inline double arma_rng_alt::randu_val () {
@@ -84,43 +83,43 @@ inline double arma_rng_alt::randn_val() {
84
83
// polar form of the Box-Muller transformation:
85
84
// http://en.wikipedia.org/wiki/Box-Muller_transformation
86
85
// http://en.wikipedia.org/wiki/Marsaglia_polar_method
87
-
86
+
88
87
double tmp1;
89
88
double tmp2;
90
89
double w;
91
-
90
+
92
91
do {
93
92
tmp1 = double (2 ) * double (::Rf_runif (0 , 1 )) - double (1 );
94
93
tmp2 = double (2 ) * double (::Rf_runif (0 , 1 )) - double (1 );
95
94
// tmp1 = double(2) * double(std::rand()) * (double(1) / double(RAND_MAX)) - double(1);
96
95
// tmp2 = double(2) * double(std::rand()) * (double(1) / double(RAND_MAX)) - double(1);
97
-
96
+
98
97
w = tmp1*tmp1 + tmp2*tmp2;
99
98
} while ( w >= double (1 ) );
100
-
99
+
101
100
return double ( tmp1 * std::sqrt ( (double (-2 ) * std::log (w)) / w) );
102
101
}
103
102
104
- template <typename eT>
103
+ template <typename eT>
105
104
inline void arma_rng_alt::randn_dual_val (eT& out1, eT& out2) {
106
105
// make sure we are internally using at least floats
107
106
typedef typename promote_type<eT,float >::result eTp;
108
-
107
+
109
108
eTp tmp1;
110
109
eTp tmp2;
111
110
eTp w;
112
-
111
+
113
112
do {
114
113
tmp1 = eTp (2 ) * eTp (::Rf_runif (0 , RAND_MAX)) * (eTp (1 ) / eTp (RAND_MAX)) - eTp (1 );
115
114
tmp2 = eTp (2 ) * eTp (::Rf_runif (0 , RAND_MAX)) * (eTp (1 ) / eTp (RAND_MAX)) - eTp (1 );
116
115
// tmp1 = eTp(2) * eTp(std::rand()) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
117
116
// tmp2 = eTp(2) * eTp(std::rand()) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
118
-
117
+
119
118
w = tmp1*tmp1 + tmp2*tmp2;
120
119
} while ( w >= eTp (1 ) );
121
-
120
+
122
121
const eTp k = std::sqrt ( (eTp (-2 ) * std::log (w)) / w);
123
-
122
+
124
123
out1 = eT (tmp1*k);
125
124
out2 = eT (tmp2*k);
126
125
}
@@ -131,14 +130,14 @@ template<typename eT>
131
130
inline void arma_rng_alt::randi_fill (eT* mem, const uword N, const int a, const int b) {
132
131
if ( (a == 0 ) && (b == RAND_MAX) ) {
133
132
for (uword i=0 ; i<N; ++i) {
134
- mem[i] = ::Rf_runif (0 , RAND_MAX);
133
+ mem[i] = static_cast <eT>( ::Rf_runif (0 , RAND_MAX)); // std::rand( );
135
134
// mem[i] = std::rand();
136
135
}
137
136
} else {
138
137
const uword length = b - a + 1 ;
139
-
138
+
140
139
const double scale = double (length) / double (RAND_MAX);
141
-
140
+
142
141
for (uword i=0 ; i<N; ++i) {
143
142
mem[i] = (std::min)( b, (int ( double (::Rf_runif (0 ,RAND_MAX)) * scale ) + a) );
144
143
// mem[i] = (std::min)( b, (int( double(std::rand()) * scale ) + a) );
0 commit comments