|
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 | 1 |
|
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/>. |
| 2 | +// This file support the legacy location and includes from the new location |
29 | 3 |
|
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 | | -classarma_rng_alt { |
42 | | -public: |
43 | | - |
44 | | -typedefunsignedint seed_type; |
45 | | - |
46 | | -inlinestaticvoidset_seed(const seed_type val); |
47 | | - |
48 | | - arma_inlinestaticintrandi_val(); |
49 | | - arma_inlinestaticdoublerandu_val(); |
50 | | -inlinestaticdoublerandn_val(); |
51 | | - |
52 | | -template<typename eT> |
53 | | -inlinestaticvoidrandn_dual_val(eT& out1, eT& out2); |
54 | | - |
55 | | -template<typename eT> |
56 | | -inlinestaticvoidrandi_fill(eT* mem,const uword N,constint a,constint b); |
57 | | - |
58 | | -inlinestaticintrandi_max_val(); |
59 | | -}; |
60 | | - |
61 | | -inlinevoidarma_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 | | -staticint 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 | | - |
74 | | -arma_inlineintarma_rng_alt::randi_val() { |
75 | | -return ::Rf_runif(0, RAND_MAX);//std::rand(); |
76 | | -} |
77 | | - |
78 | | -arma_inlinedoublearma_rng_alt::randu_val() { |
79 | | -returndouble(::Rf_runif(0,1)); |
80 | | -//return double( double(std::rand()) * ( double(1) / double(RAND_MAX) ) ); |
81 | | -} |
82 | | - |
83 | | -inlinedoublearma_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 | | -returndouble( tmp1 *std::sqrt( (double(-2) *std::log(w)) / w) ); |
102 | | -} |
103 | | - |
104 | | -template<typename eT> |
105 | | -inlinevoidarma_rng_alt::randn_dual_val(eT& out1, eT& out2) { |
106 | | -// make sure we are internally using at least floats |
107 | | -typedeftypename 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> |
131 | | -inlinevoidarma_rng_alt::randi_fill(eT* mem,const uword N,constint a,constint 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 | | -constdouble 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 | | - |
149 | | -inlineintarma_rng_alt::randi_max_val() { |
150 | | -return RAND_MAX; |
151 | | -} |
| 4 | +#include"RcppArmadillo/rng/Alt_R_RNG.h" |