Index: /library/bdm/itpp_ext.cpp
===================================================================
--- /library/bdm/itpp_ext.cpp (revision 587)
+++ /library/bdm/itpp_ext.cpp (revision 622)
@@ -18,44 +18,50 @@
 // 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,37 +74,41 @@
 #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,8 +116,10 @@
 
 //#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,21 +159,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;
 			}
@@ -173,8 +185,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;
 	}
@@ -186,19 +198,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 */
@@ -206,9 +218,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;
@@ -222,20 +234,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
@@ -244,26 +256,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;
 			}
@@ -275,12 +287,13 @@
 
 
-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,30 +301,32 @@
 	// 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,5 +336,6 @@
 #define el 0.5772156649015329
 
-double psi ( double x ) {
+double psi (double x)
+{
 	double s, ps, xa, x2;
 	int n, k;
@@ -335,67 +351,71 @@
 	};
 
-	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)
 {
Index: /applications/bdmtoolbox/tutorial/arx_mex_test.cfg
===================================================================
--- /applications/bdmtoolbox/tutorial/arx_mex_test.cfg (revision 568)
+++ /applications/bdmtoolbox/tutorial/arx_mex_test.cfg (revision 622)
@@ -1,4 +1,7 @@
 //Data generating system
-system = { type="external"; filename="arx_test.cfg"; path="system";};
+system = { type="external";
+           filename="arx_test.cfg";
+           path="system";
+         };
 
 //store results
@@ -6,18 +9,20 @@
 	type= "mexlog";
 	maxlen = 90; //
-    //dirname = "exp/ax";
+	//dirname = "exp/ax";
 };
 
 //estimation
 estimator = {
-   type = "ARX";
-	y = {type="RV"; names=("y"); };
+	type = "ARX";
+	y = {type="RV";
+	     names= ( "y" );
+	    };
 	rgr = {type="RV";
-		names = ("y","y","y","u");
-		times = [-1, -2, -3, -1];
-	};
+	       names = ( "y","y","y","u" );
+	       times = [-1, -2, -3, -1];
+	      };
 
 	//optional fields
-	dV0 = matrix(5,1,[1e-3, 1e-5, 1e-5, 1e-5, 1e-5]); //default: 1e-3 for y, 1e-5 for rgr
+	dV0 = matrix ( 5,1,[1e-3, 1e-5, 1e-5, 1e-5, 1e-5] ); //default: 1e-3 for y, 1e-5 for rgr
 	//nu0 = 8.;      //default: rgrlen + 2
 	frg = .9991;    // forgetting, default frg=1.0
@@ -25,5 +30,5 @@
 
 //experiment description
-experiment:{
+experiment: {
 	ndat = 90;
 };
Index: /applications/bdmtoolbox/tutorial/arx_test.cfg
===================================================================
--- /applications/bdmtoolbox/tutorial/arx_test.cfg (revision 357)
+++ /applications/bdmtoolbox/tutorial/arx_test.cfg (revision 622)
@@ -2,18 +2,22 @@
 system = {
 	class = "ArxDS";
-	y = {class="RV"; names=("y", "u");};
-	u = {class="RV"; names=(); };
+	y = {class="RV";
+	     names= ( "y", "u" );
+	    };
+	u = {class="RV";
+	     names=();
+	    };
 	rgr = {class="RV";
-		names = ("y","y","y","u");
-		times = [-1, -2, -3, -1];
-	};
+	       names = ( "y","y","y","u" );
+	       times = [-1, -2, -3, -1];
+	      };
 	//AR parameters
-	theta = (2,4,[0.8, -0.3, 0.4, 1.0,
-	              0.0, 0.0, 0.0, 0.0]);
+	theta = ( 2,4,[0.8, -0.3, 0.4, 1.0,
+	               0.0, 0.0, 0.0, 0.0] );
 	// offset
 	offset = [0.0, 0.0];
 	//variance
-	r = (2,2,[0.1, 0.0,
-	         0.0, 1.0] );
+	r = ( 2,2,[0.1, 0.0,
+	           0.0, 1.0] );
 	// log also theta
 	opt="L_theta";
@@ -30,9 +34,11 @@
 estimator = {
 	class = "ARX";
-	y = {class="RV"; names=("y"); };
+	y = {class="RV";
+	     names= ( "y" );
+	    };
 	rgr = {class="RV";
-		names = ("y","y","y","u");
-		times = [-1, -2, -3, -1];
-	};
+	       names = ( "y","y","y","u" );
+	       times = [-1, -2, -3, -1];
+	      };
 
 	//optional fields
@@ -43,5 +49,5 @@
 
 //experiment description
-experiment:{
+experiment: {
 	ndat = 9000;
 };
Index: /applications/bdmtoolbox/mex/simulator.cpp
===================================================================
--- /applications/bdmtoolbox/mex/simulator.cpp (revision 617)
+++ /applications/bdmtoolbox/mex/simulator.cpp (revision 622)
@@ -12,5 +12,5 @@
 \enddot
 
-Here, 
+Here,
 \li Data Source is an object (class DS) providing sequential data, \f$ [d_1, d_2, \ldots d_t] \f$.
 \li Result Logger is an object (class logger) dedicated to storing important data from the experiment.
@@ -56,23 +56,23 @@
 void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) {
 	// Check the number of inputs and output arguments
-	if(n_input<1) mexErrMsgTxt("Usage:\n"  
-		"result=simulator(system,  experiment, logger)\n"
-		"  system     = struct('class','datasource',...);  % Estimated system\n"
-		"  === optional ==="
-		"  experiment = struct('ndat',10);                 % number of data in experiment, full length of finite datasources, 10 otherwise \n"
-		"  logger     = struct('class','mexlogger');       % How to store results, default=mexlog, i.e. matlab structure\n\n"
-		"see documentation of classes datasource and mexlogger and their offsprings in BDM.");
-	
+	if ( n_input<1 ) mexErrMsgTxt ( "Usage:\n"
+		                                "result=simulator(system,  experiment, logger)\n"
+		                                "  system     = struct('class','datasource',...);  % Estimated system\n"
+		                                "  === optional ==="
+		                                "  experiment = struct('ndat',10);                 % number of data in experiment, full length of finite datasources, 10 otherwise \n"
+		                                "  logger     = struct('class','mexlogger');       % How to store results, default=mexlog, i.e. matlab structure\n\n"
+		                                "see documentation of classes datasource and mexlogger and their offsprings in BDM." );
+
 	RV::clear_all();
-	//CONFIG 
+	//CONFIG
 	UImxArray Cfg;
-	try{
-	Cfg.addGroup(input[0],"system");
-	if (n_input>2){
-		Cfg.addGroup(input[1],"experiment");
-	}
-	if (n_input>3){
-		Cfg.addGroup(input[2],"logger");
-	}/*else{
+	try {
+		Cfg.addGroup ( input[0],"system" );
+		if ( n_input>2 ) {
+			Cfg.addGroup ( input[1],"experiment" );
+		}
+		if ( n_input>3 ) {
+			Cfg.addGroup ( input[2],"logger" );
+		}/*else{
 		// define logger as mexlog
 		Setting &S=Cfg.getRoot();
@@ -85,8 +85,10 @@
 		S["logger"]["maxlen"]=maxlen;
 	}*/
-	} catch(SettingException e){it_error("error: "+string(e.getPath()));}
+	} catch ( SettingException e ) {
+		it_error ( "error: "+string ( e.getPath() ) );
+	}
 
 	//DBG
-	Cfg.writeFile("simulator.cfg");
+	Cfg.writeFile ( "simulator.cfg" );
 
 #else
@@ -100,26 +102,27 @@
 	UIFile Cfg ( fname );
 #endif
-	
-	shared_ptr<DS> Ds = UI::build<DS>( Cfg, "system" );
+
+	shared_ptr<DS> Ds = UI::build<DS> ( Cfg, "system" );
 	long Ndat=10;
-	if (Cfg.exists("experiment")) {
-		if (Cfg.lookupValue ( "experiment.ndat",Ndat )) {
-			bdm_assert(Ndat<=Ds->max_length(), "Data source has less data then required");
+	if ( Cfg.exists ( "experiment" ) ) {
+		if ( Cfg.lookupValue ( "experiment.ndat",Ndat ) ) {
+			bdm_assert ( Ndat<=Ds->max_length(), "Data source has less data then required" );
 		};
 	} else {
-		if (Ds->max_length() < std::numeric_limits< int >::max()) {
+		if ( Ds->max_length() < std::numeric_limits< int >::max() ) {
 			Ndat=Ds->max_length();
-		};// else Ndat=10;
+		}
+		;// else Ndat=10;
 	}
-	shared_ptr<logger> L = UI::build<logger>( Cfg, "logger",UI::optional);
-	if (!L) {
-		#ifdef MEX
+	shared_ptr<logger> L = UI::build<logger> ( Cfg, "logger",UI::optional );
+	if ( !L ) {
+#ifdef MEX
 		//mex logger has only from_setting constructor - we  have to fill it...
-		L=new mexlog(Ndat);
-		#else
+		L=new mexlog ( Ndat );
+#else
 		L=new stdlog();
-		#endif
+#endif
 	}
-	
+
 	Ds->log_add ( *L );
 	L->init();
@@ -130,5 +133,5 @@
 		Ds->getdata ( dt );					// read data
 		Ds->logit ( *L );
-		
+
 		L->step();
 		Ds->step();							// simulator step
@@ -139,7 +142,7 @@
 
 #ifdef MEX
-	mexlog* mL=dynamic_cast<mexlog*>(L.get());
+	mexlog* mL=dynamic_cast<mexlog*> ( L.get() );
 
-	if (mL) { // user wants output!!
+	if ( mL ) { // user wants output!!
 		if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output variables!" );
 		output[0] = mL->toCell();
Index: /applications/bdmtoolbox/mex/merger.cpp
===================================================================
--- /applications/bdmtoolbox/mex/merger.cpp (revision 569)
+++ /applications/bdmtoolbox/mex/merger.cpp (revision 622)
@@ -29,105 +29,103 @@
 
 #ifdef MEX
-void mexFunction (int n_output, mxArray *output[], int n_input, const mxArray *input[])
-{
+void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) {
 	// Check the number of inputs and output arguments
-	if (n_input != 3) mexErrMsgTxt ("Usage:\n"
-		                            "result=merger(sources, support, merger)\n"
-		                            "  sources= { struct('class','epdf'),... };  % cell of pdfs (epdfs or mpdfs) to be merged,\n"
-		                            "  support= struct(\n"
-		                            "           grid    = {[dim1_start,dim1_end], [dim2_start, dim2_end]...}  %support boundary \n"
-		                            "           nbins   = [bins_in_dim1, bins_in_dim2,...]                    %fixed \n"
-		                            "         === OR ==\n"
-		                            "           pdf     = struct('class','epdf'); % pdf to draw samples from\n"
-		                            "           nsamples= 100;                    % number of samples\n"
-		                            "           );\n"
-		                            "        If all elements are present,  (grid,nbins) is used;\n"
-		                            "  merger = struct('class','merger_*');       % object to be used for merging,\n\n"
-		                            "see documentation of classes epdf, mpdf, merger_base and their offsprings in BDM.");
+	if ( n_input != 3 ) mexErrMsgTxt ( "Usage:\n"
+		                                   "result=merger(sources, support, merger)\n"
+		                                   "  sources= { struct('class','epdf'),... };  % cell of pdfs (epdfs or mpdfs) to be merged,\n"
+		                                   "  support= struct(\n"
+		                                   "           grid    = {[dim1_start,dim1_end], [dim2_start, dim2_end]...}  %support boundary \n"
+		                                   "           nbins   = [bins_in_dim1, bins_in_dim2,...]                    %fixed \n"
+		                                   "         === OR ==\n"
+		                                   "           pdf     = struct('class','epdf'); % pdf to draw samples from\n"
+		                                   "           nsamples= 100;                    % number of samples\n"
+		                                   "           );\n"
+		                                   "        If all elements are present,  (grid,nbins) is used;\n"
+		                                   "  merger = struct('class','merger_*');       % object to be used for merging,\n\n"
+		                                   "see documentation of classes epdf, mpdf, merger_base and their offsprings in BDM." );
 	RV::clear_all();
 	// LOAD CONFIG
 	UImxArray Cfg;
-	Cfg.addList (input[0], "Sources");
-	Cfg.addGroup (input[1], "Support");
-	Cfg.addGroup (input[2], "Merger");
+	Cfg.addList ( input[0], "Sources" );
+	Cfg.addGroup ( input[1], "Support" );
+	Cfg.addGroup ( input[2], "Merger" );
 
 	//DBG
-	Cfg.writeFile ("merger.cfg");
+	Cfg.writeFile ( "merger.cfg" );
 #else
 int main() {
-	UIFile Cfg ("merger.cfg");
+	UIFile Cfg ( "merger.cfg" );
 #endif
 	// Sources
 	Array<shared_ptr<mpdf> > Sources;
 	//abuse Mer to store sources
-	Setting& _Sources = Cfg.lookup ("Sources");
+	Setting& _Sources = Cfg.lookup ( "Sources" );
 	int Slen = _Sources.getLength();
-	Sources.set_size (Slen);
-	for (int i = 0; i < Slen; i++) {
+	Sources.set_size ( Slen );
+	for ( int i = 0; i < Slen; i++ ) {
 		try {
-			shared_ptr<mpdf> mtmp = UI::build<mpdf> (_Sources, i);
-			Sources (i) = mtmp;
-		} catch (UIException) {
+			shared_ptr<mpdf> mtmp = UI::build<mpdf> ( _Sources, i );
+			Sources ( i ) = mtmp;
+		} catch ( UIException ) {
 			// it is not mpdf - see if it is epdf
 			try {
-				shared_ptr<epdf> etmp = UI::build<epdf> (_Sources, i);
-				if (etmp) {
-					Sources (i) = new mepdf (etmp); // hopefully OK
+				shared_ptr<epdf> etmp = UI::build<epdf> ( _Sources, i );
+				if ( etmp ) {
+					Sources ( i ) = new mepdf ( etmp ); // hopefully OK
 				}
-			} catch (UIException &e) {
-				it_error ("No mpdfs or epdfs found! " + string (e.what()));
-			} catch (std::exception e) {
-				it_error ("Error in UI at " + _Sources[i].getPath());
+			} catch ( UIException &e ) {
+				it_error ( "No mpdfs or epdfs found! " + string ( e.what() ) );
+			} catch ( std::exception e ) {
+				it_error ( "Error in UI at " + _Sources[i].getPath() );
 			}
-		} catch (std::exception &e) {
-			it_error ("Error in UI at " + _Sources[i].getPath());
+		} catch ( std::exception &e ) {
+			it_error ( "Error in UI at " + _Sources[i].getPath() );
 		}
 	}
 
-	shared_ptr<merger_base> Merger = UI::build<merger_base> (Cfg, "Merger");
+	shared_ptr<merger_base> Merger = UI::build<merger_base> ( Cfg, "Merger" );
 
 	// Support
 	try {
-		shared_ptr<rectangular_support> RecSup = UI::build<rectangular_support> (Cfg, "Support");
-		Merger->set_support(*RecSup);
-	}
-	catch (UIException &e) {
-		shared_ptr<discrete_support> DisSup = UI::build<discrete_support> (Cfg, "Support");
-		Merger->set_support (*DisSup);
+		shared_ptr<rectangular_support> RecSup = UI::build<rectangular_support> ( Cfg, "Support" );
+		Merger->set_support ( *RecSup );
+	} catch ( UIException &e ) {
+		shared_ptr<discrete_support> DisSup = UI::build<discrete_support> ( Cfg, "Support" );
+		Merger->set_support ( *DisSup );
 	}
 // COMPUTE RESULTS
-	Merger->set_sources (Sources); 
+	Merger->set_sources ( Sources );
 	Merger->merge();
 
 // save results
-	Array<vec> source_vals (Sources.length());
-	for (int i = 0;i < Sources.length();i++) {
+	Array<vec> source_vals ( Sources.length() );
+	for ( int i = 0;i < Sources.length();i++ ) {
 		datalink_m2e dl;
-		dl.set_connection (Sources (i)->_rv(), Sources (i)->_rvc(), Merger->_rv());
+		dl.set_connection ( Sources ( i )->_rv(), Sources ( i )->_rvc(), Merger->_rv() );
 
-		vec ll (Merger->_Smp()._samples().length());
-		for (int j = 0; j < Merger->_Smp()._samples().length(); j++) {
-			vec dt = dl.pushdown (Merger->_Smp()._samples() (j));
-			vec dtc = dl.get_cond (Merger->_Smp()._samples() (j));
-			ll (j) = Sources (i)->evallogcond (dt, dtc);
+		vec ll ( Merger->_Smp()._samples().length() );
+		for ( int j = 0; j < Merger->_Smp()._samples().length(); j++ ) {
+			vec dt = dl.pushdown ( Merger->_Smp()._samples() ( j ) );
+			vec dtc = dl.get_cond ( Merger->_Smp()._samples() ( j ) );
+			ll ( j ) = Sources ( i )->evallogcond ( dt, dtc );
 		}
 
-		vec sll = exp (ll);
+		vec sll = exp ( ll );
 
-		source_vals (i) = sll / sum (sll);
+		source_vals ( i ) = sll / sum ( sll );
 	}
 
-	merger_mix* MerMix=dynamic_cast<merger_mix*>(Merger.get());
+	merger_mix* MerMix=dynamic_cast<merger_mix*> ( Merger.get() );
 	vec mix_val;
-			
-	if (MerMix){
-		vec ll (Merger->_Smp()._samples().length());
-		for (int j = 0; j < Merger->_Smp()._samples().length(); j++) {
-			ll (j) = Merger->evallog (Merger->_Smp()._samples() (j));
+
+	if ( MerMix ) {
+		vec ll ( Merger->_Smp()._samples().length() );
+		for ( int j = 0; j < Merger->_Smp()._samples().length(); j++ ) {
+			ll ( j ) = Merger->evallog ( Merger->_Smp()._samples() ( j ) );
 		}
-	
-		vec sll = exp (ll);
-	
-		mix_val = sll / sum (sll);
+
+		vec sll = exp ( ll );
+
+		mix_val = sll / sum ( sll );
 	}
 
@@ -135,33 +133,33 @@
 	mxArray* tmp ;
 	// Save results
-	if (n_output > 0) {
-		tmp = mxCreateStructMatrix (1, 1, 0, NULL);
+	if ( n_output > 0 ) {
+		tmp = mxCreateStructMatrix ( 1, 1, 0, NULL );
 		//support
 		Array<vec> &samples = Merger->_Smp()._samples();
-		if (samples.size() > 0) {
-			mxArray* fld = mxCreateDoubleMatrix (samples (0).length(), samples.size(), mxREAL);
-			Arrayvec2mxArray (samples, fld);
-			mxReplaceFieldNM (tmp, "support", fld);
+		if ( samples.size() > 0 ) {
+			mxArray* fld = mxCreateDoubleMatrix ( samples ( 0 ).length(), samples.size(), mxREAL );
+			Arrayvec2mxArray ( samples, fld );
+			mxReplaceFieldNM ( tmp, "support", fld );
 		}
 
 		//weights
 		vec &w = Merger->_Smp()._w();
-		mxArray* fldw = mxCreateDoubleMatrix (1, w.length(), mxREAL);
-		vec2mxArray (w, fldw);
-		mxReplaceFieldNM (tmp, "weights", fldw);
+		mxArray* fldw = mxCreateDoubleMatrix ( 1, w.length(), mxREAL );
+		vec2mxArray ( w, fldw );
+		mxReplaceFieldNM ( tmp, "weights", fldw );
 
 		//mixture values
-		if (mix_val.length()>0){
-			mxArray* fldm = mxCreateDoubleMatrix (1, w.length(), mxREAL);
-			vec2mxArray (mix_val, fldm);
-			mxReplaceFieldNM (tmp, "mix", fldm);
+		if ( mix_val.length() >0 ) {
+			mxArray* fldm = mxCreateDoubleMatrix ( 1, w.length(), mxREAL );
+			vec2mxArray ( mix_val, fldm );
+			mxReplaceFieldNM ( tmp, "mix", fldm );
 		}
 		// sources
 		char srcstr[20];
-		for (int i = 0;i < Sources.length();i++) {
-			sprintf (srcstr, "source%d", i + 1);
-			mxArray* fldw = mxCreateDoubleMatrix (1, source_vals (i).length(), mxREAL);
-			vec2mxArray (source_vals (i), fldw);
-			mxReplaceFieldNM (tmp, srcstr, fldw);
+		for ( int i = 0;i < Sources.length();i++ ) {
+			sprintf ( srcstr, "source%d", i + 1 );
+			mxArray* fldw = mxCreateDoubleMatrix ( 1, source_vals ( i ).length(), mxREAL );
+			vec2mxArray ( source_vals ( i ), fldw );
+			mxReplaceFieldNM ( tmp, srcstr, fldw );
 		}
 
Index: /applications/bdmtoolbox/mex/mxstruct2config.cpp
===================================================================
--- /applications/bdmtoolbox/mex/mxstruct2config.cpp (revision 391)
+++ /applications/bdmtoolbox/mex/mxstruct2config.cpp (revision 622)
@@ -2,6 +2,6 @@
 
 void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) {
-	UImxArray F (input[0]);
-	string filename = mxArray2string(input[1]);
-	F.writeFile(filename.c_str());
+	UImxArray F ( input[0] );
+	string filename = mxArray2string ( input[1] );
+	F.writeFile ( filename.c_str() );
 }
Index: /applications/bdmtoolbox/mex/estimator.cpp
===================================================================
--- /applications/bdmtoolbox/mex/estimator.cpp (revision 611)
+++ /applications/bdmtoolbox/mex/estimator.cpp (revision 622)
@@ -9,10 +9,10 @@
 	{rank="same"; "Data Source"; "Bayesian Model"}
 	"Data Source" -> "Bayesian Model" [label="data"];
-	"Bayesian Model" -> "Result Logger" [label="estimated\n statistics"]; 
+	"Bayesian Model" -> "Result Logger" [label="estimated\n statistics"];
 	"Data Source" -> "Result Logger" [label="Simulated\n data"];
 }
 \enddot
 
-Here, 
+Here,
 \li Data Source is an object (class DS) providing sequential data, \f$ [d_1, d_2, \ldots d_t] \f$.
 \li Bayesian Model is an object (class BM) performing Bayesian filtering,
@@ -65,26 +65,26 @@
 void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) {
 	// Check the number of inputs and output arguments
-	if(n_input<2) mexErrMsgTxt("Usage:\n"  
-		"result=estimator(system, estimators, experiment, logger)\n"
-		"  system     = struct('class','datasource',...);  % Estimated system\n"
-		"  estimators = {struct('class','estimator',...),  % Estimators\n"
-		"                struct('class','estimator',...),...} \n"
-		"  === optional ==="
-		"  experiment = struct('ndat',10);                 % number of data in experiment, full length of finite datasources, 10 otherwise \n"
-		"  logger     = struct('class','mexlogger');       % How to store results, default=mexlog, i.e. matlab structure\n\n"
-		"see documentation of classes datasource, BM, and mexlogger and their offsprings in BDM.");
-	
+	if ( n_input<2 ) mexErrMsgTxt ( "Usage:\n"
+		                                "result=estimator(system, estimators, experiment, logger)\n"
+		                                "  system     = struct('class','datasource',...);  % Estimated system\n"
+		                                "  estimators = {struct('class','estimator',...),  % Estimators\n"
+		                                "                struct('class','estimator',...),...} \n"
+		                                "  === optional ==="
+		                                "  experiment = struct('ndat',10);                 % number of data in experiment, full length of finite datasources, 10 otherwise \n"
+		                                "  logger     = struct('class','mexlogger');       % How to store results, default=mexlog, i.e. matlab structure\n\n"
+		                                "see documentation of classes datasource, BM, and mexlogger and their offsprings in BDM." );
+
 	RV::clear_all();
-	//CONFIG 
+	//CONFIG
 	UImxArray Cfg;
-	try{
-	Cfg.addGroup(input[0],"system");
-	Cfg.addList(input[1],"estimators");
-	if (n_input>2){
-		Cfg.addGroup(input[2],"experiment");
-	}
-	if (n_input>3){
-		Cfg.addGroup(input[3],"logger");
-	}/*else{
+	try {
+		Cfg.addGroup ( input[0],"system" );
+		Cfg.addList ( input[1],"estimators" );
+		if ( n_input>2 ) {
+			Cfg.addGroup ( input[2],"experiment" );
+		}
+		if ( n_input>3 ) {
+			Cfg.addGroup ( input[3],"logger" );
+		}/*else{
 		// define logger as mexlog
 		Setting &S=Cfg.getRoot();
@@ -97,8 +97,10 @@
 		S["logger"]["maxlen"]=maxlen;
 	}*/
-	} catch(SettingException e){it_error("error: "+string(e.getPath()));}
+	} catch ( SettingException e ) {
+		it_error ( "error: "+string ( e.getPath() ) );
+	}
 
 	//DBG
-	Cfg.writeFile("estimator.cfg");
+	Cfg.writeFile ( "estimator.cfg" );
 
 #else
@@ -112,54 +114,56 @@
 	UIFile Cfg ( fname );
 #endif
-	
-	shared_ptr<DS> Ds = UI::build<DS>( Cfg, "system" );
-	Array<shared_ptr<BM> > Es;	UI::get(Es,Cfg, "estimators" );
+
+	shared_ptr<DS> Ds = UI::build<DS> ( Cfg, "system" );
+	Array<shared_ptr<BM> > Es;
+	UI::get ( Es,Cfg, "estimators" );
 	long Ndat=10;
-	if (Cfg.exists("experiment")) {
-		if (Cfg.lookupValue ( "experiment.ndat",Ndat )) {
-			bdm_assert(Ndat<=Ds->max_length(), "Data source has less data then required");
+	if ( Cfg.exists ( "experiment" ) ) {
+		if ( Cfg.lookupValue ( "experiment.ndat",Ndat ) ) {
+			bdm_assert ( Ndat<=Ds->max_length(), "Data source has less data then required" );
 		};
 	} else {
-		if (Ds->max_length() < inf) {
+		if ( Ds->max_length() < inf ) {
 			Ndat=Ds->max_length();
-		};// else Ndat=10;
+		}
+		;// else Ndat=10;
 	}
-	shared_ptr<logger> L = UI::build<logger>( Cfg, "logger",UI::optional);
-	if (!L) {
-		#ifdef MEX
+	shared_ptr<logger> L = UI::build<logger> ( Cfg, "logger",UI::optional );
+	if ( !L ) {
+#ifdef MEX
 		//mex logger has only from_setting constructor - we  have to fill it...
-		L=new mexlog(Ndat);
-		#else
+		L=new mexlog ( Ndat );
+#else
 		L=new stdlog();
-		#endif
+#endif
 	}
-	
+
 	Ds->log_add ( *L );
 	string Ename;
 	Setting &S=Cfg;
-	for (int i=0; i<Es.length(); i++){
-		try{
-			UI::get(Ename, S["estimators"][i], "name");
-		} catch (UIException e){
-			Ename="Est"+num2str(i);
+	for ( int i=0; i<Es.length(); i++ ) {
+		try {
+			UI::get ( Ename, S["estimators"][i], "name" );
+		} catch ( UIException e ) {
+			Ename="Est"+num2str ( i );
 		}
-		
-		Es(i)->log_add(*L,Ename); // estimate
+
+		Es ( i )->log_add ( *L,Ename ); // estimate
 	}
-		L->init();
+	L->init();
 
 	vec dt=zeros ( Ds->_drv()._dsize() );   //data variable
-	Array<datalink*> Dls(Es.length()); 
-	for (int i=0; i<Es.length(); i++){
-		Dls(i)=new datalink( Es(i)->_drv(),Ds->_drv() ); //datalink between a datasource and estimator
+	Array<datalink*> Dls ( Es.length() );
+	for ( int i=0; i<Es.length(); i++ ) {
+		Dls ( i ) =new datalink ( Es ( i )->_drv(),Ds->_drv() ); //datalink between a datasource and estimator
 	}
-	
+
 	for ( int tK=0;tK<Ndat;tK++ ) {
 		Ds->getdata ( dt );					// read data
 		Ds->logit ( *L );
-		
-		for (int i=0; i<Es.length(); i++){
-			Es(i)->bayes ( Dls(i)->pushdown ( dt ) );		// update estimates
-			Es(i)->logit (*L);
+
+		for ( int i=0; i<Es.length(); i++ ) {
+			Es ( i )->bayes ( Dls ( i )->pushdown ( dt ) );		// update estimates
+			Es ( i )->logit ( *L );
 		}
 		L->step();
@@ -171,7 +175,7 @@
 
 #ifdef MEX
-	mexlog* mL=dynamic_cast<mexlog*>(L.get());
+	mexlog* mL=dynamic_cast<mexlog*> ( L.get() );
 
-	if (mL) { // user wants output!!
+	if ( mL ) { // user wants output!!
 		if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output variables!" );
 		output[0] = mL->toCell();
