Changeset 622 for library/bdm/itpp_ext.cpp

Timestamp:

09/16/09 22:52:57 (15 years ago)

Author:

smidl

Message:

astyle

Files:

• library/bdm/itpp_ext.cpp (modified) (13 diffs)

library/bdm/itpp_ext.cpp

@@ -18 +18 @@
 // from  algebra/lapack.h
 
-extern "C" {  /* QR factorization of a general matrix A  */
-        void dgeqrf_ ( int *m, int *n, double *a, int *lda, double *tau, double *work,
-        int *lwork, int *info );
+extern "C"    /* QR factorization of a general matrix A  */
+{
+        void dgeqrf_ (int *m, int *n, double *a, int *lda, double *tau, double *work,
+                      int *lwork, int *info);
 };
 
-namespace itpp {
-Array<int> to_Arr ( const ivec &indices ) {
-        Array<int> a ( indices.size() );
-        for ( int i = 0; i < a.size(); i++ ) {
-                a ( i ) = indices ( i );
+namespace itpp
+{
+Array<int> to_Arr (const ivec &indices)
+{
+        Array<int> a (indices.size());
+        for (int i = 0; i < a.size(); i++) {
+                a (i) = indices (i);
         }
         return a;
 }
 
-ivec linspace ( int from, int to ) {
+ivec linspace (int from, int to)
+{
         int n = to - from + 1;
         int i;
-        it_assert_debug ( n > 0, "wrong linspace" );
-        ivec iv ( n );
-        for ( i = 0; i < n; i++ ) iv ( i ) = from + i;
+        it_assert_debug (n > 0, "wrong linspace");
+        ivec iv (n);
+        for (i = 0; i < n; i++) iv (i) = from + i;
         return iv;
 };
 
-void set_subvector ( vec &ov, const ivec &iv, const vec &v ) {
-        it_assert_debug ( ( iv.length() <= v.length() ),
+void set_subvector (vec &ov, const ivec &iv, const vec &v)
+{
+        it_assert_debug ( (iv.length() <= v.length()),
                           "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out "
-                          "of range of v" );
-        for ( int i = 0; i < iv.length(); i++ ) {
-                it_assert_debug ( iv ( i ) < ov.length(),
-                                  "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out "
-                                  "of range of v" );
-                ov ( iv ( i ) ) = v ( i );
-        }
-}
-
-vec get_vec ( const vec &v, const ivec &indexlist ) {
+                          "of range of v");
+        for (int i = 0; i < iv.length(); i++) {
+                it_assert_debug (iv (i) < ov.length(),
+                                 "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out "
+                                 "of range of v");
+                ov (iv (i)) = v (i);
+        }
+}
+
+vec get_vec (const vec &v, const ivec &indexlist)
+{
         int size = indexlist.size();
-        vec temp ( size );
-        for ( int i = 0; i < size; ++i ) {
-                temp ( i ) = v._data() [indexlist ( i ) ];
+        vec temp (size);
+        for (int i = 0; i < size; ++i) {
+                temp (i) = v._data() [indexlist (i) ];
         }
         return temp;
@@ -68 +74 @@
 #define R_FINITE std::isfinite
 
-bvec operator> ( const vec &t1, const vec &t2 ) {
-        it_assert_debug ( t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors" );
-        bvec temp ( t1.length() );
-        for ( int i = 0; i < t1.length(); i++ )
-                temp ( i ) = ( t1[i] > t2[i] );
+bvec operator> (const vec &t1, const vec &t2)
+{
+        it_assert_debug (t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors");
+        bvec temp (t1.length());
+        for (int i = 0; i < t1.length(); i++)
+                temp (i) = (t1[i] > t2[i]);
         return temp;
 }
 
-bvec operator< ( const vec &t1, const vec &t2 ) {
-        it_assert_debug ( t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors" );
-        bvec temp ( t1.length() );
-        for ( int i = 0; i < t1.length(); i++ )
-                temp ( i ) = ( t1[i] < t2[i] );
+bvec operator< (const vec &t1, const vec &t2)
+{
+        it_assert_debug (t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors");
+        bvec temp (t1.length());
+        for (int i = 0; i < t1.length(); i++)
+                temp (i) = (t1[i] < t2[i]);
         return temp;
 }
 
 
-bvec operator& ( const bvec &a, const bvec &b ) {
-        it_assert_debug ( b.size() == a.size(), "operator&(): Vectors of different lengths" );
-
-        bvec temp ( a.size() );
-        for ( int i = 0; i < a.size(); i++ ) {
-                temp ( i ) = a ( i ) & b ( i );
+bvec operator& (const bvec &a, const bvec &b)
+{
+        it_assert_debug (b.size() == a.size(), "operator&(): Vectors of different lengths");
+
+        bvec temp (a.size());
+        for (int i = 0; i < a.size(); i++) {
+                temp (i) = a (i) & b (i);
         }
         return temp;
 }
 
-bvec operator| ( const bvec &a, const bvec &b ) {
-        it_assert_debug ( b.size() != a.size(), "operator&(): Vectors of different lengths" );
-
-        bvec temp ( a.size() );
-        for ( int i = 0; i < a.size(); i++ ) {
-                temp ( i ) = a ( i ) | b ( i );
+bvec operator| (const bvec &a, const bvec &b)
+{
+        it_assert_debug (b.size() != a.size(), "operator&(): Vectors of different lengths");
+
+        bvec temp (a.size());
+        for (int i = 0; i < a.size(); i++) {
+                temp (i) = a (i) | b (i);
         }
         return temp;
@@ -106 +116 @@
 
 //#if 0
-Gamma_RNG::Gamma_RNG ( double a, double b ) {
-        setup ( a, b );
-}
-double  Gamma_RNG::sample() {
+Gamma_RNG::Gamma_RNG (double a, double b)
+{
+        setup (a, b);
+}
+double  Gamma_RNG::sample()
+{
         //A copy of rgamma code from the R package!!
         //
@@ -147 +159 @@
         double scale = 1.0 / beta;
 
-        if ( !R_FINITE ( a ) || !R_FINITE ( scale ) || a < 0.0 || scale <= 0.0 ) {
-                it_error ( "Gamma_RNG wrong parameters" );
-        }
-
-        if ( a < 1. ) { /* GS algorithm for parameters a < 1 */
-                if ( a == 0 )
+        if (!R_FINITE (a) || !R_FINITE (scale) || a < 0.0 || scale <= 0.0) {
+                it_error ("Gamma_RNG wrong parameters");
+        }
+
+        if (a < 1.) {   /* GS algorithm for parameters a < 1 */
+                if (a == 0)
                         return 0.;
                 e = 1.0 + exp_m1 * a;
-                for ( ;; ) {  //VS repeat
+                for (;;) {    //VS repeat
                         p = e * unif_rand();
-                        if ( p >= 1.0 ) {
-                                x = -log ( ( e - p ) / a );
-                                if ( exp_rand() >= ( 1.0 - a ) * log ( x ) )
+                        if (p >= 1.0) {
+                                x = -log ( (e - p) / a);
+                                if (exp_rand() >= (1.0 - a) * log (x))
                                         break;
                         } else {
-                                x = exp ( log ( p ) / a );
-                                if ( exp_rand() >= x )
+                                x = exp (log (p) / a);
+                                if (exp_rand() >= x)
                                         break;
                         }
@@ -173 +185 @@
 
         /* Step 1: Recalculations of s2, s, d if a has changed */
-        if ( a != aa ) {
+        if (a != aa) {
                 aa = a;
                 s2 = a - 0.5;
-                s = sqrt ( s2 );
+                s = sqrt (s2);
                 d = sqrt32 - s * 12.0;
         }
@@ -186 +198 @@
         x = s + 0.5 * t;
         ret_val = x * x;
-        if ( t >= 0.0 )
+        if (t >= 0.0)
                 return scale * ret_val;
 
         /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
         u = unif_rand();
-        if ( ( d * u ) <= ( t * t * t ) )
+        if ( (d * u) <= (t * t * t))
                 return scale * ret_val;
 
         /* Step 4: recalculations of q0, b, si, c if necessary */
 
-        if ( a != aaa ) {
+        if (a != aaa) {
                 aaa = a;
                 r = 1.0 / a;
-                q0 = ( ( ( ( ( ( q7 * r + q6 ) * r + q5 ) * r + q4 ) * r + q3 ) * r
-                         + q2 ) * r + q1 ) * r;
+                q0 = ( ( ( ( ( (q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
+                         + q2) * r + q1) * r;
 
                 /* Approximation depending on size of parameter a */
@@ -206 +218 @@
                 /* were established by numerical experiments */
 
-                if ( a <= 3.686 ) {
+                if (a <= 3.686) {
                         b = 0.463 + s + 0.178 * s2;
                         si = 1.235;
                         c = 0.195 / s - 0.079 + 0.16 * s;
-                } else if ( a <= 13.022 ) {
+                } else if (a <= 13.022) {
                         b = 1.654 + 0.0076 * s2;
                         si = 1.68 / s + 0.275;
@@ -222 +234 @@
         /* Step 5: no quotient test if x not positive */
 
-        if ( x > 0.0 ) {
+        if (x > 0.0) {
                 /* Step 6: calculation of v and quotient q */
-                v = t / ( s + s );
-                if ( fabs ( v ) <= 0.25 )
-                        q = q0 + 0.5 * t * t * ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v
-                                                     + a3 ) * v + a2 ) * v + a1 ) * v;
+                v = t / (s + s);
+                if (fabs (v) <= 0.25)
+                        q = q0 + 0.5 * t * t * ( ( ( ( ( (a7 * v + a6) * v + a5) * v + a4) * v
+                                                     + a3) * v + a2) * v + a1) * v;
                 else
-                        q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v );
+                        q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log (1.0 + v);
 
 
                 /* Step 7: quotient acceptance (q) */
-                if ( log ( 1.0 - u ) <= q )
+                if (log (1.0 - u) <= q)
                         return scale * ret_val;
         }
 
-        for ( ;; ) { //VS repeat
+        for (;;) {   //VS repeat
                 /* Step 8: e = standard exponential deviate
                  *      u =  0,1 -uniform deviate
@@ -244 +256 @@
                 u = unif_rand();
                 u = u + u - 1.0;
-                if ( u < 0.0 )
+                if (u < 0.0)
                         t = b - si * e;
                 else
                         t = b + si * e;
                 /* Step  9:  rejection if t < tau(1) = -0.71874483771719 */
-                if ( t >= -0.71874483771719 ) {
+                if (t >= -0.71874483771719) {
                         /* Step 10:      calculation of v and quotient q */
-                        v = t / ( s + s );
-                        if ( fabs ( v ) <= 0.25 )
+                        v = t / (s + s);
+                        if (fabs (v) <= 0.25)
                                 q = q0 + 0.5 * t * t *
-                                    ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v
-                                        + a2 ) * v + a1 ) * v;
+                                    ( ( ( ( ( (a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
+                                        + a2) * v + a1) * v;
                         else
-                                q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v );
+                                q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log (1.0 + v);
                         /* Step 11:      hat acceptance (h) */
                         /* (if q not positive go to step 8) */
-                        if ( q > 0.0 ) {
+                        if (q > 0.0) {
                                 // TODO: w = expm1(q);
-                                w = exp ( q ) - 1;
+                                w = exp (q) - 1;
                                 /*  ^^^^^ original code had approximation with rel.err < 2e-7 */
                                 /* if t is rejected sample again at step 8 */
-                                if ( ( c * fabs ( u ) ) <= ( w * exp ( e - 0.5 * t * t ) ) )
+                                if ( (c * fabs (u)) <= (w * exp (e - 0.5 * t * t)))
                                         break;
                         }
@@ -275 +287 @@
 
 
-bool qr ( const mat &A, mat &R ) {
+bool qr (const mat &A, mat &R)
+{
         int info;
         int m = A.rows();
         int n = A.cols();
         int lwork = n;
-        int k = std::min ( m, n );
-        vec tau ( k );
-        vec work ( lwork );
+        int k = std::min (m, n);
+        vec tau (k);
+        vec work (lwork);
 
         R = A;
@@ -288 +301 @@
         // perform workspace query for optimum lwork value
         int lwork_tmp = -1;
-        dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
-                  &info );
-        if ( info == 0 ) {
-                lwork = static_cast<int> ( work ( 0 ) );
-                work.set_size ( lwork, false );
-        }
-        dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info );
+        dgeqrf_ (&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
+                 &info);
+        if (info == 0) {
+                lwork = static_cast<int> (work (0));
+                work.set_size (lwork, false);
+        }
+        dgeqrf_ (&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
 
         // construct R
-        for ( int i = 0; i < m; i++ )
-                for ( int j = 0; j < std::min ( i, n ); j++ )
-                        R ( i, j ) = 0;
-
-        return ( info == 0 );
+        for (int i = 0; i < m; i++)
+                for (int j = 0; j < std::min (i, n); j++)
+                        R (i, j) = 0;
+
+        return (info == 0);
 }
 
 //#endif
-std::string num2str ( double d ) {
+std::string num2str (double d)
+{
         char tmp[20];//that should do
-        sprintf ( tmp, "%f", d );
-        return std::string ( tmp );
+        sprintf (tmp, "%f", d);
+        return std::string (tmp);
 };
-std::string num2str ( int i ) {
+std::string num2str (int i)
+{
         char tmp[10];//that should do
-        sprintf ( tmp, "%d", i );
-        return std::string ( tmp );
+        sprintf (tmp, "%d", i);
+        return std::string (tmp);
 };
 
@@ -321 +336 @@
 #define el 0.5772156649015329
 
-double psi ( double x ) {
+double psi (double x)
+{
         double s, ps, xa, x2;
         int n, k;
@@ -335 +351 @@
         };
 
-        xa = fabs ( x );
+        xa = fabs (x);
         s = 0.0;
-        if ( ( x == ( int ) x ) && ( x <= 0.0 ) ) {
+        if ( (x == (int) x) && (x <= 0.0)) {
                 ps = 1e308;
                 return ps;
         }
-        if ( xa == ( int ) xa ) {
+        if (xa == (int) xa) {
                 n = xa;
-                for ( k = 1; k < n; k++ ) {
+                for (k = 1; k < n; k++) {
                         s += 1.0 / k;
                 }
                 ps =  s - el;
-        } else if ( ( xa + 0.5 ) == ( ( int ) ( xa + 0.5 ) ) ) {
+        } else if ( (xa + 0.5) == ( (int) (xa + 0.5))) {
                 n = xa - 0.5;
-                for ( k = 1; k <= n; k++ ) {
-                        s += 1.0 / ( 2.0 * k - 1.0 );
+                for (k = 1; k <= n; k++) {
+                        s += 1.0 / (2.0 * k - 1.0);
                 }
                 ps = 2.0 * s - el - 1.386294361119891;
         } else {
-                if ( xa < 10.0 ) {
-                        n = 10 - ( int ) xa;
-                        for ( k = 0; k < n; k++ ) {
-                                s += 1.0 / ( xa + k );
+                if (xa < 10.0) {
+                        n = 10 - (int) xa;
+                        for (k = 0; k < n; k++) {
+                                s += 1.0 / (xa + k);
                         }
                         xa += n;
                 }
-                x2 = 1.0 / ( xa * xa );
-                ps = log ( xa ) - 0.5 / xa + x2 * ( ( ( ( ( ( ( a[7] * x2 + a[6] ) * x2 + a[5] ) * x2 +
-                                                        a[4] ) * x2 + a[3] ) * x2 + a[2] ) * x2 + a[1] ) * x2 + a[0] );
+                x2 = 1.0 / (xa * xa);
+                ps = log (xa) - 0.5 / xa + x2 * ( ( ( ( ( ( (a[7] * x2 + a[6]) * x2 + a[5]) * x2 +
+                                                        a[4]) * x2 + a[3]) * x2 + a[2]) * x2 + a[1]) * x2 + a[0]);
                 ps -= s;
         }
-        if ( x < 0.0 )
-                ps = ps - M_PI * std::cos ( M_PI * x ) / std::sin ( M_PI * x ) - 1.0 / x;
+        if (x < 0.0)
+                ps = ps - M_PI * std::cos (M_PI * x) / std::sin (M_PI * x) - 1.0 / x;
         return ps;
 }
 
-void triu(mat &A){
-        for(int i=1;i<A.rows();i++) { // row cycle
-                for (int j=0; j<i; j++) {A(i,j)=0;}
-        }
-}
-
-class RandunStorage{
-        const int A;
-        const int M;
-        static double seed;
-        static int counter;
+void triu (mat &A)
+{
+        for (int i = 1;i < A.rows();i++) { // row cycle
+                for (int j = 0; j < i; j++) {A (i, j) = 0;}
+        }
+}
+
+class RandunStorage
+{
+                const int A;
+                const int M;
+                static double seed;
+                static int counter;
         public:
-                RandunStorage(): A(16807), M(2147483647) {};
-                void set_seed(double seed0){seed=seed0;}
-                double get(){
+                RandunStorage() : A (16807), M (2147483647) {};
+                void set_seed (double seed0) {seed = seed0;}
+                double get() {
                         long long tmp = A * seed;
                         tmp = tmp % M;
                         seed = tmp;
-                        counter++; 
-                        return seed/M;}
+                        counter++;
+                        return seed / M;
+                }
 };
 static RandunStorage randun_global_storage;
-double RandunStorage::seed=1111111;
-int RandunStorage::counter=0;
-double randun(){return randun_global_storage.get();};
-vec randun(int n){vec res(n); for(int i=0;i<n;i++){res(i)=randun();}; return res;};
-mat randun(int n, int m){mat res(n,m); for(int i=0;i<n*m;i++){res(i)=randun();}; return res;};
+double RandunStorage::seed = 1111111;
+int RandunStorage::counter = 0;
+double randun() {return randun_global_storage.get();};
+vec randun (int n) {vec res (n); for (int i = 0;i < n;i++) {res (i) = randun();}; return res;};
+mat randun (int n, int m) {mat res (n, m); for (int i = 0;i < n*m;i++) {res (i) = randun();}; return res;};
+
 ivec unique (const ivec &in)
 {