1 | // |
---|
2 | // C++ Implementation: itpp_ext |
---|
3 | // |
---|
4 | // Description: |
---|
5 | // |
---|
6 | // |
---|
7 | // Author: smidl <smidl@utia.cas.cz>, (C) 2008 |
---|
8 | // |
---|
9 | // Copyright: See COPYING file that comes with this distribution |
---|
10 | // |
---|
11 | // |
---|
12 | |
---|
13 | #include "itpp_ext.h" |
---|
14 | |
---|
15 | // from algebra/lapack.h |
---|
16 | |
---|
17 | extern "C" { /* QR factorization of a general matrix A */ |
---|
18 | void dgeqrf_(int *m, int *n, double *a, int *lda, double *tau, double *work, |
---|
19 | int *lwork, int *info); |
---|
20 | }; |
---|
21 | |
---|
22 | namespace itpp { |
---|
23 | Array<int> to_Arr(const ivec &indices) |
---|
24 | { |
---|
25 | Array<int> a(indices.size()); |
---|
26 | for (int i = 0; i < a.size(); i++) { |
---|
27 | a(i) = indices(i); |
---|
28 | } |
---|
29 | return a; |
---|
30 | } |
---|
31 | |
---|
32 | ivec linspace ( int from, int to ) { |
---|
33 | int n=to-from+1; |
---|
34 | int i; |
---|
35 | it_assert_debug ( n>0,"wrong linspace" ); |
---|
36 | ivec iv ( n ); for (i=0;i<n;i++) iv(i)=from+i; |
---|
37 | return iv; |
---|
38 | }; |
---|
39 | |
---|
40 | //Gamma |
---|
41 | |
---|
42 | Gamma_RNG::Gamma_RNG(double a, double b) |
---|
43 | { |
---|
44 | setup(a,b); |
---|
45 | } |
---|
46 | |
---|
47 | #define log std::log |
---|
48 | #define exp std::exp |
---|
49 | #define sqrt std::sqrt |
---|
50 | #define R_FINITE std::isfinite |
---|
51 | |
---|
52 | double Gamma_RNG::sample() |
---|
53 | { |
---|
54 | //A copy of rgamma code from the R package!! |
---|
55 | // |
---|
56 | |
---|
57 | /* Constants : */ |
---|
58 | const static double sqrt32 = 5.656854; |
---|
59 | const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */ |
---|
60 | |
---|
61 | /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k)) |
---|
62 | * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k) |
---|
63 | * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k) |
---|
64 | */ |
---|
65 | const static double q1 = 0.04166669; |
---|
66 | const static double q2 = 0.02083148; |
---|
67 | const static double q3 = 0.00801191; |
---|
68 | const static double q4 = 0.00144121; |
---|
69 | const static double q5 = -7.388e-5; |
---|
70 | const static double q6 = 2.4511e-4; |
---|
71 | const static double q7 = 2.424e-4; |
---|
72 | |
---|
73 | const static double a1 = 0.3333333; |
---|
74 | const static double a2 = -0.250003; |
---|
75 | const static double a3 = 0.2000062; |
---|
76 | const static double a4 = -0.1662921; |
---|
77 | const static double a5 = 0.1423657; |
---|
78 | const static double a6 = -0.1367177; |
---|
79 | const static double a7 = 0.1233795; |
---|
80 | |
---|
81 | /* State variables [FIXME for threading!] :*/ |
---|
82 | static double aa = 0.; |
---|
83 | static double aaa = 0.; |
---|
84 | static double s, s2, d; /* no. 1 (step 1) */ |
---|
85 | static double q0, b, si, c;/* no. 2 (step 4) */ |
---|
86 | |
---|
87 | double e, p, q, r, t, u, v, w, x, ret_val; |
---|
88 | double a=alpha; |
---|
89 | double scale=1.0/beta; |
---|
90 | |
---|
91 | if (!R_FINITE(a) || !R_FINITE(scale) || a < 0.0 || scale <= 0.0) |
---|
92 | it_error("Gamma_RNG wrong parameters"); |
---|
93 | |
---|
94 | if (a < 1.) { /* GS algorithm for parameters a < 1 */ |
---|
95 | if(a == 0) |
---|
96 | return 0.; |
---|
97 | e = 1.0 + exp_m1 * a; |
---|
98 | for(;;) { //VS repeat |
---|
99 | p = e * unif_rand(); |
---|
100 | if (p >= 1.0) { |
---|
101 | x = -log((e - p) / a); |
---|
102 | if (exp_rand() >= (1.0 - a) * log(x)) |
---|
103 | break; |
---|
104 | } else { |
---|
105 | x = exp(log(p) / a); |
---|
106 | if (exp_rand() >= x) |
---|
107 | break; |
---|
108 | } |
---|
109 | } |
---|
110 | return scale * x; |
---|
111 | } |
---|
112 | |
---|
113 | /* --- a >= 1 : GD algorithm --- */ |
---|
114 | |
---|
115 | /* Step 1: Recalculations of s2, s, d if a has changed */ |
---|
116 | if (a != aa) { |
---|
117 | aa = a; |
---|
118 | s2 = a - 0.5; |
---|
119 | s = sqrt(s2); |
---|
120 | d = sqrt32 - s * 12.0; |
---|
121 | } |
---|
122 | /* Step 2: t = standard normal deviate, |
---|
123 | x = (s,1/2) -normal deviate. */ |
---|
124 | |
---|
125 | /* immediate acceptance (i) */ |
---|
126 | t = norm_rand(); |
---|
127 | x = s + 0.5 * t; |
---|
128 | ret_val = x * x; |
---|
129 | if (t >= 0.0) |
---|
130 | return scale * ret_val; |
---|
131 | |
---|
132 | /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */ |
---|
133 | u = unif_rand(); |
---|
134 | if ((d * u) <= (t * t * t)) |
---|
135 | return scale * ret_val; |
---|
136 | |
---|
137 | /* Step 4: recalculations of q0, b, si, c if necessary */ |
---|
138 | |
---|
139 | if (a != aaa) { |
---|
140 | aaa = a; |
---|
141 | r = 1.0 / a; |
---|
142 | q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r |
---|
143 | + q2) * r + q1) * r; |
---|
144 | |
---|
145 | /* Approximation depending on size of parameter a */ |
---|
146 | /* The constants in the expressions for b, si and c */ |
---|
147 | /* were established by numerical experiments */ |
---|
148 | |
---|
149 | if (a <= 3.686) { |
---|
150 | b = 0.463 + s + 0.178 * s2; |
---|
151 | si = 1.235; |
---|
152 | c = 0.195 / s - 0.079 + 0.16 * s; |
---|
153 | } else if (a <= 13.022) { |
---|
154 | b = 1.654 + 0.0076 * s2; |
---|
155 | si = 1.68 / s + 0.275; |
---|
156 | c = 0.062 / s + 0.024; |
---|
157 | } else { |
---|
158 | b = 1.77; |
---|
159 | si = 0.75; |
---|
160 | c = 0.1515 / s; |
---|
161 | } |
---|
162 | } |
---|
163 | /* Step 5: no quotient test if x not positive */ |
---|
164 | |
---|
165 | if (x > 0.0) { |
---|
166 | /* Step 6: calculation of v and quotient q */ |
---|
167 | v = t / (s + s); |
---|
168 | if (fabs(v) <= 0.25) |
---|
169 | q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v |
---|
170 | + a3) * v + a2) * v + a1) * v; |
---|
171 | else |
---|
172 | q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v); |
---|
173 | |
---|
174 | |
---|
175 | /* Step 7: quotient acceptance (q) */ |
---|
176 | if (log(1.0 - u) <= q) |
---|
177 | return scale * ret_val; |
---|
178 | } |
---|
179 | |
---|
180 | for(;;) { //VS repeat |
---|
181 | /* Step 8: e = standard exponential deviate |
---|
182 | * u = 0,1 -uniform deviate |
---|
183 | * t = (b,si)-double exponential (laplace) sample */ |
---|
184 | e = exp_rand(); |
---|
185 | u = unif_rand(); |
---|
186 | u = u + u - 1.0; |
---|
187 | if (u < 0.0) |
---|
188 | t = b - si * e; |
---|
189 | else |
---|
190 | t = b + si * e; |
---|
191 | /* Step 9: rejection if t < tau(1) = -0.71874483771719 */ |
---|
192 | if (t >= -0.71874483771719) { |
---|
193 | /* Step 10: calculation of v and quotient q */ |
---|
194 | v = t / (s + s); |
---|
195 | if (fabs(v) <= 0.25) |
---|
196 | q = q0 + 0.5 * t * t * |
---|
197 | ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v |
---|
198 | + a2) * v + a1) * v; |
---|
199 | else |
---|
200 | q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v); |
---|
201 | /* Step 11: hat acceptance (h) */ |
---|
202 | /* (if q not positive go to step 8) */ |
---|
203 | if (q > 0.0) { |
---|
204 | // TODO: w = expm1(q); |
---|
205 | w = exp(q)-1; |
---|
206 | /* ^^^^^ original code had approximation with rel.err < 2e-7 */ |
---|
207 | /* if t is rejected sample again at step 8 */ |
---|
208 | if ((c * fabs(u)) <= (w * exp(e - 0.5 * t * t))) |
---|
209 | break; |
---|
210 | } |
---|
211 | } |
---|
212 | } /* repeat .. until `t' is accepted */ |
---|
213 | x = s + 0.5 * t; |
---|
214 | return scale * x * x; |
---|
215 | } |
---|
216 | |
---|
217 | |
---|
218 | bool qr(const mat &A, mat &R) |
---|
219 | { |
---|
220 | int info; |
---|
221 | int m = A.rows(); |
---|
222 | int n = A.cols(); |
---|
223 | int lwork = n; |
---|
224 | int k = std::min(m, n); |
---|
225 | vec tau(k); |
---|
226 | vec work(lwork); |
---|
227 | |
---|
228 | R = A; |
---|
229 | |
---|
230 | // perform workspace query for optimum lwork value |
---|
231 | int lwork_tmp = -1; |
---|
232 | dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp, |
---|
233 | &info); |
---|
234 | if (info == 0) { |
---|
235 | lwork = static_cast<int>(work(0)); |
---|
236 | work.set_size(lwork, false); |
---|
237 | } |
---|
238 | dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info); |
---|
239 | |
---|
240 | // construct R |
---|
241 | for (int i=0; i<m; i++) |
---|
242 | for (int j=0; j<std::min(i,n); j++) |
---|
243 | R(i,j) = 0; |
---|
244 | |
---|
245 | return (info==0); |
---|
246 | } |
---|
247 | |
---|
248 | } |
---|