Index: library/bdm/itpp_ext.cpp
===================================================================
--- library/bdm/itpp_ext.cpp (revision 723)
+++ library/bdm/itpp_ext.cpp (revision 737)
@@ -18,52 +18,46 @@
 // 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);
-};
-
-namespace itpp
-{
-	vec empty_vec = vec(0);
-	
-Array<int> to_Arr (const ivec &indices)
-{
-	Array<int> a (indices.size());
-	for (int i = 0; i < a.size(); i++) {
-		a (i) = indices (i);
+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 {
+vec empty_vec = vec ( 0 );
+
+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;
@@ -76,41 +70,37 @@
 #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]);
-	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]);
-	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);
-	}
-	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 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] );
+	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 );
+	}
+	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 );
 	}
 	return temp;
@@ -118,25 +108,23 @@
 
 //! poor man's operator vec(bvec) - copied for svn version of itpp
-ivec get_from_bvec(const ivec &v, const bvec &binlist){
-  int size = binlist.size();
-  it_assert_debug(v.size() == size, "Vec<>::operator()(bvec &): "
-                  "Wrong size of binlist vector");
-  ivec temp(size);
-  int j = 0;
-  for (int i = 0; i < size; ++i)
-    if (binlist(i) == bin(1))
-      temp(j++) = v(i);
-  temp.set_size(j, true);
-  return temp;
+ivec get_from_bvec ( const ivec &v, const bvec &binlist ) {
+	int size = binlist.size();
+	it_assert_debug ( v.size() == size, "Vec<>::operator()(bvec &): "
+	                  "Wrong size of binlist vector" );
+	ivec temp ( size );
+	int j = 0;
+	for ( int i = 0; i < size; ++i )
+		if ( binlist ( i ) == bin ( 1 ) )
+			temp ( j++ ) = v ( i );
+	temp.set_size ( j, true );
+	return temp;
 }
 
 
 //#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!!
 	//
@@ -176,21 +164,21 @@
 	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;
 			}
@@ -202,8 +190,8 @@
 
 	/* 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;
 	}
@@ -215,19 +203,19 @@
 	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 */
@@ -235,9 +223,9 @@
 		/* 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;
@@ -251,20 +239,20 @@
 	/* 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
@@ -273,26 +261,26 @@
 		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;
 			}
@@ -304,13 +292,12 @@
 
 
-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;
@@ -318,32 +305,30 @@
 	// 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);
-};
-std::string num2str (int i)
-{
+	sprintf ( tmp, "%f", d );
+	return std::string ( tmp );
+};
+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 );
 };
 
@@ -353,6 +338,5 @@
 #define el 0.5772156649015329
 
-double psi (double x)
-{
+double psi ( double x ) {
 	double s, ps, xa, x2;
 	int n, k;
@@ -368,110 +352,124 @@
 	};
 
-	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;}
+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;
+		}
 	}
 }
 
 //! Storage of randun() internals
-class RandunStorage
-{
-		const int A;
-		const int M;
-		static double seed;
-		static int counter;
-	public:
-		RandunStorage() : A (16807), M (2147483647) {};
-		//!set seed of the randun() generator
-		void set_seed (double seed0) {seed = seed0;}
-		//! generate randun() sample
-		double get() {
-			long long tmp = A * seed;
-			tmp = tmp % M;
-			seed = tmp;
-			counter++;
-			return seed / M;
-		}
+class RandunStorage {
+	const int A;
+	const int M;
+	static double seed;
+	static int counter;
+public:
+	RandunStorage() : A ( 16807 ), M ( 2147483647 ) {};
+	//!set seed of the randun() generator
+	void set_seed ( double seed0 ) {
+		seed = seed0;
+	}
+	//! generate randun() sample
+	double get() {
+		long long tmp = A * seed;
+		tmp = tmp % M;
+		seed = tmp;
+		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;};
-
-ivec unique (const ivec &in)
-{
-	ivec uniq (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 ) {
+	ivec uniq ( 0 );
 	int j = 0;
 	bool found = false;
-	for (int i = 0;i < in.length(); i++) {
+	for ( int i = 0; i < in.length(); i++ ) {
 		found = false;
 		j = 0;
-		while ( (!found) && (j < uniq.length())) {
-			if (in (i) == uniq (j)) found = true;
+		while ( ( !found ) && ( j < uniq.length() ) ) {
+			if ( in ( i ) == uniq ( j ) ) found = true;
 			j++;
 		}
-		if (!found) uniq = concat (uniq, in (i));
+		if ( !found ) uniq = concat ( uniq, in ( i ) );
 	}
 	return uniq;
 }
 
-ivec unique_complement (const ivec &in, const ivec &base)
-{
+ivec unique_complement ( const ivec &in, const ivec &base ) {
 	// almost a copy of unique
-	ivec uniq (0);
+	ivec uniq ( 0 );
 	int j = 0;
 	bool found = false;
-	for (int i = 0;i < in.length(); i++) {
+	for ( int i = 0; i < in.length(); i++ ) {
 		found = false;
 		j = 0;
-		while ( (!found) && (j < uniq.length())) {
-			if (in (i) == uniq (j)) found = true;
+		while ( ( !found ) && ( j < uniq.length() ) ) {
+			if ( in ( i ) == uniq ( j ) ) found = true;
 			j++;
 		}
-		j=0;
-		while ( (!found) && (j < base.length())) {
-			if (in (i) == base (j)) found = true;
+		j = 0;
+		while ( ( !found ) && ( j < base.length() ) ) {
+			if ( in ( i ) == base ( j ) ) found = true;
 			j++;
 		}
-		if (!found) uniq = concat (uniq, in (i));
+		if ( !found ) uniq = concat ( uniq, in ( i ) );
 	}
 	return uniq;
