Changeset 737 for library/bdm/itpp_ext.cpp
- Timestamp:
- 11/25/09 12:14:38 (15 years ago)
- Files:
-
- 1 modified
Legend:
- Unmodified
- Added
- Removed
-
library/bdm/itpp_ext.cpp
r723 r737 18 18 // from algebra/lapack.h 19 19 20 extern "C" /* QR factorization of a general matrix A */ 21 { 22 void dgeqrf_ (int *m, int *n, double *a, int *lda, double *tau, double *work, 23 int *lwork, int *info); 24 }; 25 26 namespace itpp 27 { 28 vec empty_vec = vec(0); 29 30 Array<int> to_Arr (const ivec &indices) 31 { 32 Array<int> a (indices.size()); 33 for (int i = 0; i < a.size(); i++) { 34 a (i) = indices (i); 20 extern "C" { /* QR factorization of a general matrix A */ 21 void dgeqrf_ ( int *m, int *n, double *a, int *lda, double *tau, double *work, 22 int *lwork, int *info ); 23 }; 24 25 namespace itpp { 26 vec empty_vec = vec ( 0 ); 27 28 Array<int> to_Arr ( const ivec &indices ) { 29 Array<int> a ( indices.size() ); 30 for ( int i = 0; i < a.size(); i++ ) { 31 a ( i ) = indices ( i ); 35 32 } 36 33 return a; 37 34 } 38 35 39 ivec linspace (int from, int to) 40 { 36 ivec linspace ( int from, int to ) { 41 37 int n = to - from + 1; 42 38 int i; 43 it_assert_debug ( n > 0, "wrong linspace");44 ivec iv ( n);45 for ( i = 0; i < n; i++) iv (i) = from + i;39 it_assert_debug ( n > 0, "wrong linspace" ); 40 ivec iv ( n ); 41 for ( i = 0; i < n; i++ ) iv ( i ) = from + i; 46 42 return iv; 47 43 }; 48 44 49 void set_subvector (vec &ov, const ivec &iv, const vec &v) 50 { 51 it_assert_debug ( (iv.length() <= v.length()), 45 void set_subvector ( vec &ov, const ivec &iv, const vec &v ) { 46 it_assert_debug ( ( iv.length() <= v.length() ), 52 47 "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out " 53 "of range of v"); 54 for (int i = 0; i < iv.length(); i++) { 55 it_assert_debug (iv (i) < ov.length(), 56 "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out " 57 "of range of v"); 58 ov (iv (i)) = v (i); 59 } 60 } 61 62 vec get_vec (const vec &v, const ivec &indexlist) 63 { 48 "of range of v" ); 49 for ( int i = 0; i < iv.length(); i++ ) { 50 it_assert_debug ( iv ( i ) < ov.length(), 51 "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out " 52 "of range of v" ); 53 ov ( iv ( i ) ) = v ( i ); 54 } 55 } 56 57 vec get_vec ( const vec &v, const ivec &indexlist ) { 64 58 int size = indexlist.size(); 65 vec temp ( size);66 for ( int i = 0; i < size; ++i) {67 temp ( i) = v._data() [indexlist (i) ];59 vec temp ( size ); 60 for ( int i = 0; i < size; ++i ) { 61 temp ( i ) = v._data() [indexlist ( i ) ]; 68 62 } 69 63 return temp; … … 76 70 #define R_FINITE std::isfinite 77 71 78 bvec operator> (const vec &t1, const vec &t2) 79 { 80 it_assert_debug (t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors"); 81 bvec temp (t1.length()); 82 for (int i = 0; i < t1.length(); i++) 83 temp (i) = (t1[i] > t2[i]); 84 return temp; 85 } 86 87 bvec operator< (const vec &t1, const vec &t2) 88 { 89 it_assert_debug (t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors"); 90 bvec temp (t1.length()); 91 for (int i = 0; i < t1.length(); i++) 92 temp (i) = (t1[i] < t2[i]); 93 return temp; 94 } 95 96 97 bvec operator& (const bvec &a, const bvec &b) 98 { 99 it_assert_debug (b.size() == a.size(), "operator&(): Vectors of different lengths"); 100 101 bvec temp (a.size()); 102 for (int i = 0; i < a.size(); i++) { 103 temp (i) = a (i) & b (i); 104 } 105 return temp; 106 } 107 108 bvec operator| (const bvec &a, const bvec &b) 109 { 110 it_assert_debug (b.size() == a.size(), "operator|(): Vectors of different lengths"); 111 112 bvec temp (a.size()); 113 for (int i = 0; i < a.size(); i++) { 114 temp (i) = a (i) | b (i); 72 bvec operator> ( const vec &t1, const vec &t2 ) { 73 it_assert_debug ( t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors" ); 74 bvec temp ( t1.length() ); 75 for ( int i = 0; i < t1.length(); i++ ) 76 temp ( i ) = ( t1[i] > t2[i] ); 77 return temp; 78 } 79 80 bvec operator< ( const vec &t1, const vec &t2 ) { 81 it_assert_debug ( t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors" ); 82 bvec temp ( t1.length() ); 83 for ( int i = 0; i < t1.length(); i++ ) 84 temp ( i ) = ( t1[i] < t2[i] ); 85 return temp; 86 } 87 88 89 bvec operator& ( const bvec &a, const bvec &b ) { 90 it_assert_debug ( b.size() == a.size(), "operator&(): Vectors of different lengths" ); 91 92 bvec temp ( a.size() ); 93 for ( int i = 0; i < a.size(); i++ ) { 94 temp ( i ) = a ( i ) & b ( i ); 95 } 96 return temp; 97 } 98 99 bvec operator| ( const bvec &a, const bvec &b ) { 100 it_assert_debug ( b.size() == a.size(), "operator|(): Vectors of different lengths" ); 101 102 bvec temp ( a.size() ); 103 for ( int i = 0; i < a.size(); i++ ) { 104 temp ( i ) = a ( i ) | b ( i ); 115 105 } 116 106 return temp; … … 118 108 119 109 //! poor man's operator vec(bvec) - copied for svn version of itpp 120 ivec get_from_bvec (const ivec &v, const bvec &binlist){121 122 it_assert_debug(v.size() == size, "Vec<>::operator()(bvec &): "123 "Wrong size of binlist vector");124 ivec temp(size);125 126 for (int i = 0; i < size; ++i)127 if (binlist(i) == bin(1))128 temp(j++) = v(i);129 temp.set_size(j, true);130 110 ivec get_from_bvec ( const ivec &v, const bvec &binlist ) { 111 int size = binlist.size(); 112 it_assert_debug ( v.size() == size, "Vec<>::operator()(bvec &): " 113 "Wrong size of binlist vector" ); 114 ivec temp ( size ); 115 int j = 0; 116 for ( int i = 0; i < size; ++i ) 117 if ( binlist ( i ) == bin ( 1 ) ) 118 temp ( j++ ) = v ( i ); 119 temp.set_size ( j, true ); 120 return temp; 131 121 } 132 122 133 123 134 124 //#if 0 135 Gamma_RNG::Gamma_RNG (double a, double b) 136 { 137 setup (a, b); 138 } 139 double Gamma_RNG::sample() 140 { 125 Gamma_RNG::Gamma_RNG ( double a, double b ) { 126 setup ( a, b ); 127 } 128 double Gamma_RNG::sample() { 141 129 //A copy of rgamma code from the R package!! 142 130 // … … 176 164 double scale = 1.0 / beta; 177 165 178 if ( !R_FINITE (a) || !R_FINITE (scale) || a < 0.0 || scale <= 0.0) {179 it_error ( "Gamma_RNG wrong parameters");180 } 181 182 if ( a < 1.) {/* GS algorithm for parameters a < 1 */183 if ( a == 0)166 if ( !R_FINITE ( a ) || !R_FINITE ( scale ) || a < 0.0 || scale <= 0.0 ) { 167 it_error ( "Gamma_RNG wrong parameters" ); 168 } 169 170 if ( a < 1. ) { /* GS algorithm for parameters a < 1 */ 171 if ( a == 0 ) 184 172 return 0.; 185 173 e = 1.0 + exp_m1 * a; 186 for ( ;;) {//VS repeat174 for ( ;; ) { //VS repeat 187 175 p = e * unif_rand(); 188 if ( p >= 1.0) {189 x = -log ( ( e - p) / a);190 if ( exp_rand() >= (1.0 - a) * log (x))176 if ( p >= 1.0 ) { 177 x = -log ( ( e - p ) / a ); 178 if ( exp_rand() >= ( 1.0 - a ) * log ( x ) ) 191 179 break; 192 180 } else { 193 x = exp ( log (p) / a);194 if ( exp_rand() >= x)181 x = exp ( log ( p ) / a ); 182 if ( exp_rand() >= x ) 195 183 break; 196 184 } … … 202 190 203 191 /* Step 1: Recalculations of s2, s, d if a has changed */ 204 if ( a != aa) {192 if ( a != aa ) { 205 193 aa = a; 206 194 s2 = a - 0.5; 207 s = sqrt ( s2);195 s = sqrt ( s2 ); 208 196 d = sqrt32 - s * 12.0; 209 197 } … … 215 203 x = s + 0.5 * t; 216 204 ret_val = x * x; 217 if ( t >= 0.0)205 if ( t >= 0.0 ) 218 206 return scale * ret_val; 219 207 220 208 /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */ 221 209 u = unif_rand(); 222 if ( ( d * u) <= (t * t * t))210 if ( ( d * u ) <= ( t * t * t ) ) 223 211 return scale * ret_val; 224 212 225 213 /* Step 4: recalculations of q0, b, si, c if necessary */ 226 214 227 if ( a != aaa) {215 if ( a != aaa ) { 228 216 aaa = a; 229 217 r = 1.0 / a; 230 q0 = ( ( ( ( ( ( q7 * r + q6) * r + q5) * r + q4) * r + q3) * r231 + q2 ) * r + q1) * r;218 q0 = ( ( ( ( ( ( q7 * r + q6 ) * r + q5 ) * r + q4 ) * r + q3 ) * r 219 + q2 ) * r + q1 ) * r; 232 220 233 221 /* Approximation depending on size of parameter a */ … … 235 223 /* were established by numerical experiments */ 236 224 237 if ( a <= 3.686) {225 if ( a <= 3.686 ) { 238 226 b = 0.463 + s + 0.178 * s2; 239 227 si = 1.235; 240 228 c = 0.195 / s - 0.079 + 0.16 * s; 241 } else if ( a <= 13.022) {229 } else if ( a <= 13.022 ) { 242 230 b = 1.654 + 0.0076 * s2; 243 231 si = 1.68 / s + 0.275; … … 251 239 /* Step 5: no quotient test if x not positive */ 252 240 253 if ( x > 0.0) {241 if ( x > 0.0 ) { 254 242 /* Step 6: calculation of v and quotient q */ 255 v = t / ( s + s);256 if ( fabs (v) <= 0.25)257 q = q0 + 0.5 * t * t * ( ( ( ( ( ( a7 * v + a6) * v + a5) * v + a4) * v258 + a3 ) * v + a2) * v + a1) * v;243 v = t / ( s + s ); 244 if ( fabs ( v ) <= 0.25 ) 245 q = q0 + 0.5 * t * t * ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v 246 + a3 ) * v + a2 ) * v + a1 ) * v; 259 247 else 260 q = q0 - s * t + 0.25 * t * t + ( s2 + s2) * log (1.0 + v);248 q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v ); 261 249 262 250 263 251 /* Step 7: quotient acceptance (q) */ 264 if ( log (1.0 - u) <= q)252 if ( log ( 1.0 - u ) <= q ) 265 253 return scale * ret_val; 266 254 } 267 255 268 for ( ;;) {//VS repeat256 for ( ;; ) { //VS repeat 269 257 /* Step 8: e = standard exponential deviate 270 258 * u = 0,1 -uniform deviate … … 273 261 u = unif_rand(); 274 262 u = u + u - 1.0; 275 if ( u < 0.0)263 if ( u < 0.0 ) 276 264 t = b - si * e; 277 265 else 278 266 t = b + si * e; 279 267 /* Step 9: rejection if t < tau(1) = -0.71874483771719 */ 280 if ( t >= -0.71874483771719) {268 if ( t >= -0.71874483771719 ) { 281 269 /* Step 10: calculation of v and quotient q */ 282 v = t / ( s + s);283 if ( fabs (v) <= 0.25)270 v = t / ( s + s ); 271 if ( fabs ( v ) <= 0.25 ) 284 272 q = q0 + 0.5 * t * t * 285 ( ( ( ( ( ( a7 * v + a6) * v + a5) * v + a4) * v + a3) * v286 + a2 ) * v + a1) * v;273 ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v 274 + a2 ) * v + a1 ) * v; 287 275 else 288 q = q0 - s * t + 0.25 * t * t + ( s2 + s2) * log (1.0 + v);276 q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v ); 289 277 /* Step 11: hat acceptance (h) */ 290 278 /* (if q not positive go to step 8) */ 291 if ( q > 0.0) {279 if ( q > 0.0 ) { 292 280 // TODO: w = expm1(q); 293 w = exp ( q) - 1;281 w = exp ( q ) - 1; 294 282 /* ^^^^^ original code had approximation with rel.err < 2e-7 */ 295 283 /* if t is rejected sample again at step 8 */ 296 if ( ( c * fabs (u)) <= (w * exp (e - 0.5 * t * t)))284 if ( ( c * fabs ( u ) ) <= ( w * exp ( e - 0.5 * t * t ) ) ) 297 285 break; 298 286 } … … 304 292 305 293 306 bool qr (const mat &A, mat &R) 307 { 294 bool qr ( const mat &A, mat &R ) { 308 295 int info; 309 296 int m = A.rows(); 310 297 int n = A.cols(); 311 298 int lwork = n; 312 int k = std::min ( m, n);313 vec tau ( k);314 vec work ( lwork);299 int k = std::min ( m, n ); 300 vec tau ( k ); 301 vec work ( lwork ); 315 302 316 303 R = A; … … 318 305 // perform workspace query for optimum lwork value 319 306 int lwork_tmp = -1; 320 dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,321 &info);322 if ( info == 0) {323 lwork = static_cast<int> ( work (0));324 work.set_size ( lwork, false);325 } 326 dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);307 dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp, 308 &info ); 309 if ( info == 0 ) { 310 lwork = static_cast<int> ( work ( 0 ) ); 311 work.set_size ( lwork, false ); 312 } 313 dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info ); 327 314 328 315 // construct R 329 for ( int i = 0; i < m; i++)330 for ( int j = 0; j < std::min (i, n); j++)331 R ( i, j) = 0;332 333 return ( info == 0);316 for ( int i = 0; i < m; i++ ) 317 for ( int j = 0; j < std::min ( i, n ); j++ ) 318 R ( i, j ) = 0; 319 320 return ( info == 0 ); 334 321 } 335 322 336 323 //#endif 337 std::string num2str (double d) 338 { 324 std::string num2str ( double d ) { 339 325 char tmp[20];//that should do 340 sprintf (tmp, "%f", d); 341 return std::string (tmp); 342 }; 343 std::string num2str (int i) 344 { 326 sprintf ( tmp, "%f", d ); 327 return std::string ( tmp ); 328 }; 329 std::string num2str ( int i ) { 345 330 char tmp[10];//that should do 346 sprintf ( tmp, "%d", i);347 return std::string ( tmp);331 sprintf ( tmp, "%d", i ); 332 return std::string ( tmp ); 348 333 }; 349 334 … … 353 338 #define el 0.5772156649015329 354 339 355 double psi (double x) 356 { 340 double psi ( double x ) { 357 341 double s, ps, xa, x2; 358 342 int n, k; … … 368 352 }; 369 353 370 xa = fabs ( x);354 xa = fabs ( x ); 371 355 s = 0.0; 372 if ( ( x == (int) x) && (x <= 0.0)) {356 if ( ( x == ( int ) x ) && ( x <= 0.0 ) ) { 373 357 ps = 1e308; 374 358 return ps; 375 359 } 376 if ( xa == (int) xa) {360 if ( xa == ( int ) xa ) { 377 361 n = xa; 378 for ( k = 1; k < n; k++) {362 for ( k = 1; k < n; k++ ) { 379 363 s += 1.0 / k; 380 364 } 381 365 ps = s - el; 382 } else if ( ( xa + 0.5) == ( (int) (xa + 0.5))) {366 } else if ( ( xa + 0.5 ) == ( ( int ) ( xa + 0.5 ) ) ) { 383 367 n = xa - 0.5; 384 for ( k = 1; k <= n; k++) {385 s += 1.0 / ( 2.0 * k - 1.0);368 for ( k = 1; k <= n; k++ ) { 369 s += 1.0 / ( 2.0 * k - 1.0 ); 386 370 } 387 371 ps = 2.0 * s - el - 1.386294361119891; 388 372 } else { 389 if ( xa < 10.0) {390 n = 10 - ( int) xa;391 for ( k = 0; k < n; k++) {392 s += 1.0 / ( xa + k);373 if ( xa < 10.0 ) { 374 n = 10 - ( int ) xa; 375 for ( k = 0; k < n; k++ ) { 376 s += 1.0 / ( xa + k ); 393 377 } 394 378 xa += n; 395 379 } 396 x2 = 1.0 / ( xa * xa);397 ps = log ( xa) - 0.5 / xa + x2 * ( ( ( ( ( ( (a[7] * x2 + a[6]) * x2 + a[5]) * x2 +398 a[4] ) * x2 + a[3]) * x2 + a[2]) * x2 + a[1]) * x2 + a[0]);380 x2 = 1.0 / ( xa * xa ); 381 ps = log ( xa ) - 0.5 / xa + x2 * ( ( ( ( ( ( ( a[7] * x2 + a[6] ) * x2 + a[5] ) * x2 + 382 a[4] ) * x2 + a[3] ) * x2 + a[2] ) * x2 + a[1] ) * x2 + a[0] ); 399 383 ps -= s; 400 384 } 401 if ( x < 0.0)402 ps = ps - M_PI * std::cos ( M_PI * x) / std::sin (M_PI * x) - 1.0 / x;385 if ( x < 0.0 ) 386 ps = ps - M_PI * std::cos ( M_PI * x ) / std::sin ( M_PI * x ) - 1.0 / x; 403 387 return ps; 404 388 } 405 389 406 void triu (mat &A) 407 { 408 for (int i = 1;i < A.rows();i++) { // row cycle 409 for (int j = 0; j < i; j++) {A (i, j) = 0;} 390 void triu ( mat &A ) { 391 for ( int i = 1; i < A.rows(); i++ ) { // row cycle 392 for ( int j = 0; j < i; j++ ) { 393 A ( i, j ) = 0; 394 } 410 395 } 411 396 } 412 397 413 398 //! Storage of randun() internals 414 class RandunStorage 415 { 416 const int A; 417 const int M; 418 static double seed; 419 static int counter; 420 public: 421 RandunStorage() : A (16807), M (2147483647) {}; 422 //!set seed of the randun() generator 423 void set_seed (double seed0) {seed = seed0;} 424 //! generate randun() sample 425 double get() { 426 long long tmp = A * seed; 427 tmp = tmp % M; 428 seed = tmp; 429 counter++; 430 return seed / M; 431 } 399 class RandunStorage { 400 const int A; 401 const int M; 402 static double seed; 403 static int counter; 404 public: 405 RandunStorage() : A ( 16807 ), M ( 2147483647 ) {}; 406 //!set seed of the randun() generator 407 void set_seed ( double seed0 ) { 408 seed = seed0; 409 } 410 //! generate randun() sample 411 double get() { 412 long long tmp = A * seed; 413 tmp = tmp % M; 414 seed = tmp; 415 counter++; 416 return seed / M; 417 } 432 418 }; 433 419 static RandunStorage randun_global_storage; 434 420 double RandunStorage::seed = 1111111; 435 421 int RandunStorage::counter = 0; 436 double randun() {return randun_global_storage.get();}; 437 vec randun (int n) {vec res (n); for (int i = 0;i < n;i++) {res (i) = randun();}; return res;}; 438 mat randun (int n, int m) {mat res (n, m); for (int i = 0;i < n*m;i++) {res (i) = randun();}; return res;}; 439 440 ivec unique (const ivec &in) 441 { 442 ivec uniq (0); 422 double randun() { 423 return randun_global_storage.get(); 424 }; 425 vec randun ( int n ) { 426 vec res ( n ); 427 for ( int i = 0; i < n; i++ ) { 428 res ( i ) = randun(); 429 }; 430 return res; 431 }; 432 mat randun ( int n, int m ) { 433 mat res ( n, m ); 434 for ( int i = 0; i < n*m; i++ ) { 435 res ( i ) = randun(); 436 }; 437 return res; 438 }; 439 440 ivec unique ( const ivec &in ) { 441 ivec uniq ( 0 ); 443 442 int j = 0; 444 443 bool found = false; 445 for ( int i = 0;i < in.length(); i++) {444 for ( int i = 0; i < in.length(); i++ ) { 446 445 found = false; 447 446 j = 0; 448 while ( ( !found) && (j < uniq.length())) {449 if ( in (i) == uniq (j)) found = true;447 while ( ( !found ) && ( j < uniq.length() ) ) { 448 if ( in ( i ) == uniq ( j ) ) found = true; 450 449 j++; 451 450 } 452 if ( !found) uniq = concat (uniq, in (i));451 if ( !found ) uniq = concat ( uniq, in ( i ) ); 453 452 } 454 453 return uniq; 455 454 } 456 455 457 ivec unique_complement (const ivec &in, const ivec &base) 458 { 456 ivec unique_complement ( const ivec &in, const ivec &base ) { 459 457 // almost a copy of unique 460 ivec uniq ( 0);458 ivec uniq ( 0 ); 461 459 int j = 0; 462 460 bool found = false; 463 for ( int i = 0;i < in.length(); i++) {461 for ( int i = 0; i < in.length(); i++ ) { 464 462 found = false; 465 463 j = 0; 466 while ( ( !found) && (j < uniq.length())) {467 if ( in (i) == uniq (j)) found = true;464 while ( ( !found ) && ( j < uniq.length() ) ) { 465 if ( in ( i ) == uniq ( j ) ) found = true; 468 466 j++; 469 467 } 470 j =0;471 while ( ( !found) && (j < base.length())) {472 if ( in (i) == base (j)) found = true;468 j = 0; 469 while ( ( !found ) && ( j < base.length() ) ) { 470 if ( in ( i ) == base ( j ) ) found = true; 473 471 j++; 474 472 } 475 if ( !found) uniq = concat (uniq, in (i));473 if ( !found ) uniq = concat ( uniq, in ( i ) ); 476 474 } 477 475 return uniq;