Index: library/bdm/stat/merger.cpp
===================================================================
--- library/bdm/stat/merger.cpp (revision 461)
+++ library/bdm/stat/merger.cpp (revision 477)
@@ -3,252 +3,237 @@
 #include "../estim/arx.h"
 
-namespace bdm
-{
-
-	merger_base::merger_base(const Array<mpdf*> &S, bool own) {
-	    DBG = false;
-	    dbg_file = NULL;
-	    set_sources(S, own);
-	}
-
-	vec merger_base::merge_points ( mat &lW ) {
-		int nu=lW.rows();
-		vec result;
-		ivec indW;
-		bool infexist;
-		switch ( METHOD ) {
-			case ARITHMETIC:
-				result= log ( sum ( exp ( lW ) ) ); //ugly!
-				break;
-			case GEOMETRIC:
-				result= sum ( lW ) /nu;
-				break;
-			case LOGNORMAL:
-				vec sumlW=sum ( lW ) ;
-				indW=find((sumlW<inf) & (sumlW>-inf));
-				infexist=(indW.size()<lW.cols());
-				vec mu;
-				vec lam;
-				if (infexist){
-					mu = sumlW(indW) /nu; //mean of logs
-					// 
-					mat validlW=lW.get_cols(indW);
-					lam = sum ( pow ( validlW-outer_product ( ones ( validlW.rows() ),mu ),2 ) );
-				}
-				else {
-					mu = sum ( lW ) /nu; //mean of logs
-					lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) );
-				}
-				//
-				double coef=0.0;
-				vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere
-				switch ( nu ) {
-					case 2:
-						coef= ( 1-0.5*sqrt ( ( 4.0*beta-3.0 ) /beta ) );
-						result =coef*sq2bl + mu ;
-						break;
-					// case 4: == can be done similar to case 2 - is it worth it???
-					default: // see accompanying document merge_lognorm_derivation.lyx
-						coef = sqrt(1-(nu+1)/(2*beta*nu));
-						result =log(besselk((nu-3)/2, sq2bl*coef))-log(besselk((nu-3)/2, sq2bl)) + mu;
-						break;
-				}
-				break;
-		}
-		if (infexist){
-			vec tmp =-inf*ones(lW.cols());
-			set_subvector(tmp, indW, result);
-			return tmp;
-		}
-		else {return result;}
-	}
-
-	void merger_mix::merge ( )
-	{
-		Array<vec> &Smp = eSmp._samples(); //aux
-		vec &w = eSmp._w(); //aux
-
-		mat Smp_ex =ones ( dim +1,Npoints ); // Extended samples for the ARX model - the last row is ones
-		for ( int i=0;i<Npoints;i++ ) {	set_col_part ( Smp_ex,i,Smp ( i ) );}
-
-		if ( DBG )	*dbg_file << Name ( "Smp_0" ) << Smp_ex;
-
-		// Stuff for merging
-		vec lw_src ( Npoints );		// weights of the ith source
-		vec lw_mix ( Npoints );		// weights of the approximating mixture
-		vec lw ( Npoints );			// tmp
-		mat lW=zeros ( Nsources,Npoints ); // array of weights of all sources
-		vec vec0 ( 0 );
-
-		//initialize importance weights
-		lw_mix = 1.0; // assuming uniform grid density -- otherwise 
-
-		// Initial component in the mixture model
-		mat V0=1e-8*eye ( dim +1 );
-		ARX A0;	A0.set_statistics ( dim, V0 ); //initial guess of Mix:
-
-		Mix.init ( &A0, Smp_ex, Ncoms );
-		//Preserve initial mixture for repetitive estimation via flattening
-		MixEF Mix_init ( Mix );
-
-		// ============= MAIN LOOP ==================
-		bool converged=false;
-		int niter = 0;
-		char dbg_str[100];
-
-		emix* Mpred=Mix.epredictor ( );
-		vec Mix_pdf ( Npoints );
-		while ( !converged )
-		{
-			//Re-estimate Mix
-			//Re-Initialize Mixture model
-			Mix.flatten ( &Mix_init );
-			Mix.bayesB ( Smp_ex, w*Npoints );
-			delete Mpred;
-			Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!!
-			Mpred->set_rv ( rv ); //the predictor predicts rv of this merger
-
-			// This will be active only later in iterations!!!
-			if (  1./sum_sqr ( w ) <effss_coef*Npoints ) 
-			{
-				// Generate new samples
-				eSmp.set_samples ( Mpred );
-				for ( int i=0;i<Npoints;i++ )
-				{
-					//////////// !!!!!!!!!!!!!
-					//if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; }
-					set_col_part ( Smp_ex,i,Smp ( i ) );
-					//Importance of the mixture
-					//lw_mix ( i ) =Mix.logpred (Smp_ex.get_col(i) );
-					lw_mix ( i ) = Mpred->evallog ( Smp ( i ) );
-				}
-				if ( DBG )
-				{
-					cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;
-					cout << Mix._e()->mean() <<endl;
-					cout << sum ( Smp_ex,2 ) /Npoints <<endl;
-					cout << Smp_ex*Smp_ex.T() /Npoints << endl;
-				}
-			}
-			if ( DBG )
-			{
-				sprintf ( dbg_str,"Mpred_mean%d",niter );
-				*dbg_file << Name ( dbg_str ) << Mpred->mean();
-				sprintf ( dbg_str,"Mpred_var%d",niter );
-				*dbg_file << Name ( dbg_str ) << Mpred->variance();
-
-
-				sprintf ( dbg_str,"Mpdf%d",niter );
-				for ( int i=0;i<Npoints;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );}
-				*dbg_file << Name ( dbg_str ) << Mix_pdf;
-
-				sprintf ( dbg_str,"Smp%d",niter );
-				*dbg_file << Name ( dbg_str ) << Smp_ex;
-
-			}
-			//Importace weighting
-			for ( int i=0;i<mpdfs.length();i++ )
-			{
-				lw_src=0.0;
-				//======== Same RVs ===========
-				//Split according to dependency in rvs
-				if ( mpdfs ( i )->dimension() ==dim )
-				{
-					// no need for conditioning or marginalization
-					lw_src = mpdfs ( i )->e()->evallog_m ( Smp  );
-				}
-				else
-				{
-					// compute likelihood of marginal on the conditional variable
-					if ( mpdfs ( i )->dimensionc() >0 )
-					{
-						// Make marginal on rvc_i
-						epdf* tmp_marg = Mpred->marginal ( mpdfs ( i )->_rvc() );
-						//compute vector of lw_src
-						for ( int k=0;k<Npoints;k++ )
-						{
-							// Here val of tmp_marg = cond of mpdfs(i) ==> calling dls->get_cond
-							lw_src ( k ) += tmp_marg->evallog ( dls ( i )->get_cond ( Smp ( k ) ) );
-						}
-						delete tmp_marg;
+namespace bdm {
+
+merger_base::merger_base ( const Array<mpdf*> &S, bool own ) {
+	DBG = false;
+	dbg_file = NULL;
+	set_sources ( S, own );
+}
+
+vec merger_base::merge_points ( mat &lW ) {
+	int nu = lW.rows();
+	vec result;
+	ivec indW;
+	bool infexist;
+	switch ( METHOD ) {
+	case ARITHMETIC:
+		result = log ( sum ( exp ( lW ) ) ); //ugly!
+		break;
+	case GEOMETRIC:
+		result = sum ( lW ) / nu;
+		break;
+	case LOGNORMAL:
+		vec sumlW = sum ( lW ) ;
+		indW = find ( ( sumlW < inf ) & ( sumlW > -inf ) );
+		infexist = ( indW.size() < lW.cols() );
+		vec mu;
+		vec lam;
+		if ( infexist ) {
+			mu = sumlW ( indW ) / nu; //mean of logs
+			//
+			mat validlW = lW.get_cols ( indW );
+			lam = sum ( pow ( validlW - outer_product ( ones ( validlW.rows() ), mu ), 2 ) );
+		} else {
+			mu = sum ( lW ) / nu; //mean of logs
+			lam = sum ( pow ( lW - outer_product ( ones ( lW.rows() ), mu ), 2 ) );
+		}
+		//
+		double coef = 0.0;
+		vec sq2bl = sqrt ( 2 * beta * lam ); //this term is everywhere
+		switch ( nu ) {
+		case 2:
+			coef = ( 1 - 0.5 * sqrt ( ( 4.0 * beta - 3.0 ) / beta ) );
+			result = coef * sq2bl + mu ;
+			break;
+			// case 4: == can be done similar to case 2 - is it worth it???
+		default: // see accompanying document merge_lognorm_derivation.lyx
+			coef = sqrt ( 1 - ( nu + 1 ) / ( 2 * beta * nu ) );
+			result = log ( besselk ( ( nu - 3 ) / 2, sq2bl * coef ) ) - log ( besselk ( ( nu - 3 ) / 2, sq2bl ) ) + mu;
+			break;
+		}
+		break;
+	}
+	if ( infexist ) {
+		vec tmp = -inf * ones ( lW.cols() );
+		set_subvector ( tmp, indW, result );
+		return tmp;
+	} else {
+		return result;
+	}
+}
+
+void merger_mix::merge ( ) {
+	Array<vec> &Smp = eSmp._samples(); //aux
+	vec &w = eSmp._w(); //aux
+
+	mat Smp_ex = ones ( dim + 1, Npoints ); // Extended samples for the ARX model - the last row is ones
+	for ( int i = 0; i < Npoints; i++ ) {
+		set_col_part ( Smp_ex, i, Smp ( i ) );
+	}
+
+	if ( DBG )	*dbg_file << Name ( "Smp_0" ) << Smp_ex;
+
+	// Stuff for merging
+	vec lw_src ( Npoints );		// weights of the ith source
+	vec lw_mix ( Npoints );		// weights of the approximating mixture
+	vec lw ( Npoints );			// tmp
+	mat lW = zeros ( Nsources, Npoints ); // array of weights of all sources
+	vec vec0 ( 0 );
+
+	//initialize importance weights
+	lw_mix = 1.0; // assuming uniform grid density -- otherwise
+
+	// Initial component in the mixture model
+	mat V0 = 1e-8 * eye ( dim + 1 );
+	ARX A0;
+	A0.set_statistics ( dim, V0 ); //initial guess of Mix:
+
+	Mix.init ( &A0, Smp_ex, Ncoms );
+	//Preserve initial mixture for repetitive estimation via flattening
+	MixEF Mix_init ( Mix );
+
+	// ============= MAIN LOOP ==================
+	bool converged = false;
+	int niter = 0;
+	char dbg_str[100];
+
+	emix* Mpred = Mix.epredictor ( );
+	vec Mix_pdf ( Npoints );
+	while ( !converged ) {
+		//Re-estimate Mix
+		//Re-Initialize Mixture model
+		Mix.flatten ( &Mix_init );
+		Mix.bayesB ( Smp_ex, w*Npoints );
+		delete Mpred;
+		Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!!
+		Mpred->set_rv ( rv ); //the predictor predicts rv of this merger
+
+		// This will be active only later in iterations!!!
+		if ( 1. / sum_sqr ( w ) < effss_coef*Npoints ) {
+			// Generate new samples
+			eSmp.set_samples ( Mpred );
+			for ( int i = 0; i < Npoints; i++ ) {
+				//////////// !!!!!!!!!!!!!
+				//if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; }
+				set_col_part ( Smp_ex, i, Smp ( i ) );
+				//Importance of the mixture
+				//lw_mix ( i ) =Mix.logpred (Smp_ex.get_col(i) );
+				lw_mix ( i ) = Mpred->evallog ( Smp ( i ) );
+			}
+			if ( DBG ) {
+				cout << "Resampling =" << 1. / sum_sqr ( w ) << endl;
+				cout << Mix._e()->mean() << endl;
+				cout << sum ( Smp_ex, 2 ) / Npoints << endl;
+				cout << Smp_ex*Smp_ex.T() / Npoints << endl;
+			}
+		}
+		if ( DBG ) {
+			sprintf ( dbg_str, "Mpred_mean%d", niter );
+			*dbg_file << Name ( dbg_str ) << Mpred->mean();
+			sprintf ( dbg_str, "Mpred_var%d", niter );
+			*dbg_file << Name ( dbg_str ) << Mpred->variance();
+
+
+			sprintf ( dbg_str, "Mpdf%d", niter );
+			for ( int i = 0; i < Npoints; i++ ) {
+				Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );
+			}
+			*dbg_file << Name ( dbg_str ) << Mix_pdf;
+
+			sprintf ( dbg_str, "Smp%d", niter );
+			*dbg_file << Name ( dbg_str ) << Smp_ex;
+
+		}
+		//Importace weighting
+		for ( int i = 0; i < mpdfs.length(); i++ ) {
+			lw_src = 0.0;
+			//======== Same RVs ===========
+			//Split according to dependency in rvs
+			if ( mpdfs ( i )->dimension() == dim ) {
+				// no need for conditioning or marginalization
+				lw_src = mpdfs ( i )->e()->evallog_m ( Smp );
+			} else {
+				// compute likelihood of marginal on the conditional variable
+				if ( mpdfs ( i )->dimensionc() > 0 ) {
+					// Make marginal on rvc_i
+					epdf* tmp_marg = Mpred->marginal ( mpdfs ( i )->_rvc() );
+					//compute vector of lw_src
+					for ( int k = 0; k < Npoints; k++ ) {
+						// Here val of tmp_marg = cond of mpdfs(i) ==> calling dls->get_cond
+						lw_src ( k ) += tmp_marg->evallog ( dls ( i )->get_cond ( Smp ( k ) ) );
+					}
+					delete tmp_marg;
 
 // 					sprintf ( str,"marg%d",niter );
 // 					*dbg << Name ( str ) << lw_src;
 
+				}
+				// Compute likelihood of the missing variable
+				if ( dim > ( mpdfs ( i )->dimension() + mpdfs ( i )->dimensionc() ) ) {
+					///////////////
+					// There are variales unknown to mpdfs(i) : rvzs
+					mpdf* tmp_cond = Mpred->condition ( rvzs ( i ) );
+					// Compute likelihood
+					vec lw_dbg = lw_src;
+					for ( int k = 0; k < Npoints; k++ ) {
+						lw_src ( k ) += log (
+						                    tmp_cond->evallogcond (
+						                        zdls ( i )->pushdown ( Smp ( k ) ),
+						                        zdls ( i )->get_cond ( Smp ( k ) ) ) );
+						if ( !std::isfinite ( lw_src ( k ) ) ) {
+							lw_src ( k ) = -1e16;
+							cout << "!";
+						}
 					}
-					// Compute likelihood of the missing variable
-					if ( dim > ( mpdfs ( i )->dimension() + mpdfs ( i )->dimensionc() ) )
-					{
-						///////////////
-						// There are variales unknown to mpdfs(i) : rvzs
-						mpdf* tmp_cond = Mpred->condition ( rvzs ( i ) );
-						// Compute likelihood
-						vec lw_dbg=lw_src;
-						for ( int k= 0; k<Npoints; k++ )
-						{
-							lw_src ( k ) += log (
-							                    tmp_cond->evallogcond (
-							                        zdls ( i )->pushdown ( Smp ( k ) ),
-							                        zdls ( i )->get_cond ( Smp ( k ) ) ) );
-							if ( !std::isfinite ( lw_src ( k ) ) )
-							{
-								lw_src ( k ) = -1e16; cout << "!";
-							}
-						}
-						delete tmp_cond;
-					}
-					// Compute likelihood of the partial source
-					for ( int k= 0; k<Npoints; k++ )
-					{
-						mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp ( k ) ) );
-						lw_src ( k ) += mpdfs ( i )->e()->evallog ( dls ( i )->pushdown ( Smp ( k ) ) );
-					}
-
+					delete tmp_cond;
 				}
-		//			it_assert_debug(std::isfinite(sum(lw_src)),"bad");
-				lW.set_row ( i, lw_src ); // do not divide by mix
-			}
-			lw = merger_base::merge_points ( lW ); //merge
-
-			//Importance weighting
-			lw -=  lw_mix; // hoping that it is not numerically sensitive...
-			w = exp ( lw-max ( lw ) );
-
-			//renormalize
-			double sumw =sum ( w );
-			if ( std::isfinite ( sumw ) )
-			{
-				w = w/sumw;
-			}
-			else
-			{
-				it_file itf ( "merg_err.it" );
-				itf << Name ( "w" ) << w;
-			}
-
-			if ( DBG )
-			{
-				sprintf ( dbg_str,"lW%d",niter );
-				*dbg_file << Name ( dbg_str ) << lW;
-				sprintf ( dbg_str,"w%d",niter );
-				*dbg_file << Name ( dbg_str ) << w;
-				sprintf ( dbg_str,"lw_m%d",niter );
-				*dbg_file << Name ( dbg_str ) << lw_mix;
-			}
-			// ==== stopping rule ===
-			niter++;
-			converged = ( niter>stop_niter );
-		}
-		delete Mpred;
+				// Compute likelihood of the partial source
+				for ( int k = 0; k < Npoints; k++ ) {
+					mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp ( k ) ) );
+					lw_src ( k ) += mpdfs ( i )->e()->evallog ( dls ( i )->pushdown ( Smp ( k ) ) );
+				}
+
+			}
+			//			it_assert_debug(std::isfinite(sum(lw_src)),"bad");
+			lW.set_row ( i, lw_src ); // do not divide by mix
+		}
+		lw = merger_base::merge_points ( lW ); //merge
+
+		//Importance weighting
+		lw -=  lw_mix; // hoping that it is not numerically sensitive...
+		w = exp ( lw - max ( lw ) );
+
+		//renormalize
+		double sumw = sum ( w );
+		if ( std::isfinite ( sumw ) ) {
+			w = w / sumw;
+		} else {
+			it_file itf ( "merg_err.it" );
+			itf << Name ( "w" ) << w;
+		}
+
+		if ( DBG ) {
+			sprintf ( dbg_str, "lW%d", niter );
+			*dbg_file << Name ( dbg_str ) << lW;
+			sprintf ( dbg_str, "w%d", niter );
+			*dbg_file << Name ( dbg_str ) << w;
+			sprintf ( dbg_str, "lw_m%d", niter );
+			*dbg_file << Name ( dbg_str ) << lw_mix;
+		}
+		// ==== stopping rule ===
+		niter++;
+		converged = ( niter > stop_niter );
+	}
+	delete Mpred;
 //		cout << endl;
 
-	}
-
-	// DEFAULTS FOR MERGER_BASE 
-	const MERGER_METHOD merger_base::DFLT_METHOD=LOGNORMAL;
-	const double merger_base::DFLT_beta=1.2;
-	// DEFAULTS FOR MERGER_MIX
-	const int merger_mix::DFLT_Ncoms=10;
-	const double merger_mix::DFLT_effss_coef=0.5;
-
-}
+}
+
+// DEFAULTS FOR MERGER_BASE
+const MERGER_METHOD merger_base::DFLT_METHOD = LOGNORMAL;
+const double merger_base::DFLT_beta = 1.2;
+// DEFAULTS FOR MERGER_MIX
+const int merger_mix::DFLT_Ncoms = 10;
+const double merger_mix::DFLT_effss_coef = 0.5;
+
+}
Index: library/bdm/stat/emix.cpp
===================================================================
--- library/bdm/stat/emix.cpp (revision 464)
+++ library/bdm/stat/emix.cpp (revision 477)
@@ -1,29 +1,30 @@
 #include "emix.h"
 
-namespace bdm{
+namespace bdm {
 
 void emix::set_parameters ( const vec &w0, const Array<epdf*> &Coms0, bool copy ) {
-	w = w0/sum ( w0 );
-	dim = Coms0(0)->dimension();
-	bool isnamed = Coms0(0)->isnamed();
+	w = w0 / sum ( w0 );
+	dim = Coms0 ( 0 )->dimension();
+	bool isnamed = Coms0 ( 0 )->isnamed();
 	int i;
 	RV tmp_rv;
-	if (isnamed) tmp_rv=Coms0(0)->_rv();
-	
-	for ( i=0;i<w.length();i++ ) {
-		it_assert_debug ( dim== ( Coms0 ( i )->dimension() ),"Component sizes do not match!" );	
-		it_assert_debug ( !isnamed || tmp_rv.equal( Coms0 ( i )->_rv() ),"Component RVs do not match!" );
+	if ( isnamed ) tmp_rv = Coms0 ( 0 )->_rv();
+
+	for ( i = 0; i < w.length(); i++ ) {
+		it_assert_debug ( dim == ( Coms0 ( i )->dimension() ), "Component sizes do not match!" );
+		it_assert_debug ( !isnamed || tmp_rv.equal ( Coms0 ( i )->_rv() ), "Component RVs do not match!" );
 	}
 	if ( copy ) {
-		Coms.set_length(Coms0.length());
-		for ( i=0;i<w.length();i++ ) {it_error("Not imp...");
-			*Coms ( i ) =*Coms0 ( i );}
-		destroyComs=true;
-	}
-	else {
+		Coms.set_length ( Coms0.length() );
+		for ( i = 0; i < w.length(); i++ ) {
+			it_error ( "Not imp..." );
+			*Coms ( i ) = *Coms0 ( i );
+		}
+		destroyComs = true;
+	} else {
 		Coms = Coms0;
-		destroyComs=false;
-	}
-	if (isnamed) epdf::set_rv(tmp_rv); //coms aer already OK, no need for set_rv
+		destroyComs = false;
+	}
+	if ( isnamed ) epdf::set_rv ( tmp_rv ); //coms aer already OK, no need for set_rv
 }
 
@@ -32,45 +33,50 @@
 	vec cumDist = cumsum ( w );
 	double u0;
-	#pragma omp critical
+#pragma omp critical
 	u0 = UniRNG.sample();
 
-	int i=0;
-	while ( ( cumDist ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
+	int i = 0;
+	while ( ( cumDist ( i ) < u0 ) && ( i < ( w.length() - 1 ) ) ) {
+		i++;
+	}
 
 	return Coms ( i )->sample();
 }
 
-emix* emix::marginal(const RV &rv) const{
-	it_assert_debug(isnamed(), "rvs are not assigned");
-			
-	Array<epdf*> Cn(Coms.length());
-	for(int i=0;i<Coms.length();i++){Cn(i)=Coms(i)->marginal(rv);}
+emix* emix::marginal ( const RV &rv ) const {
+	it_assert_debug ( isnamed(), "rvs are not assigned" );
+
+	Array<epdf*> Cn ( Coms.length() );
+	for ( int i = 0; i < Coms.length(); i++ ) {
+		Cn ( i ) = Coms ( i )->marginal ( rv );
+	}
 	emix* tmp = new emix();
-	tmp->set_parameters(w,Cn,false);
+	tmp->set_parameters ( w, Cn, false );
 	tmp->ownComs();
 	return tmp;
 }
 
-mratio* emix::condition(const RV &rv) const{
-	it_assert_debug(isnamed(), "rvs are not assigned");
-	return new mratio(this,rv);
+mratio* emix::condition ( const RV &rv ) const {
+	it_assert_debug ( isnamed(), "rvs are not assigned" );
+	return new mratio ( this, rv );
 };
 
 void egiwmix::set_parameters ( const vec &w0, const Array<egiw*> &Coms0, bool copy ) {
-	w = w0/sum ( w0 );
-	dim = Coms0(0)->dimension();
+	w = w0 / sum ( w0 );
+	dim = Coms0 ( 0 )->dimension();
 	int i;
-	for ( i=0;i<w.length();i++ ) {
-		it_assert_debug ( dim== ( Coms0 ( i )->dimension() ),"Component sizes do not match!" );
+	for ( i = 0; i < w.length(); i++ ) {
+		it_assert_debug ( dim == ( Coms0 ( i )->dimension() ), "Component sizes do not match!" );
 	}
 	if ( copy ) {
-		Coms.set_length(Coms0.length());
-		for ( i=0;i<w.length();i++ ) {it_error("Not imp...");
-			*Coms ( i ) =*Coms0 ( i );}
-		destroyComs=true;
-	}
-	else {
+		Coms.set_length ( Coms0.length() );
+		for ( i = 0; i < w.length(); i++ ) {
+			it_error ( "Not imp..." );
+			*Coms ( i ) = *Coms0 ( i );
+		}
+		destroyComs = true;
+	} else {
 		Coms = Coms0;
-		destroyComs=false;
+		destroyComs = false;
 	}
 }
@@ -80,9 +86,11 @@
 	vec cumDist = cumsum ( w );
 	double u0;
-	#pragma omp critical
+#pragma omp critical
 	u0 = UniRNG.sample();
 
-	int i=0;
-	while ( ( cumDist ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
+	int i = 0;
+	while ( ( cumDist ( i ) < u0 ) && ( i < ( w.length() - 1 ) ) ) {
+		i++;
+	}
 
 	return Coms ( i )->sample();
@@ -90,6 +98,9 @@
 
 vec egiwmix::mean() const {
-	int i; vec mu = zeros ( dim );
-	for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
+	int i;
+	vec mu = zeros ( dim );
+	for ( i = 0; i < w.length(); i++ ) {
+		mu += w ( i ) * Coms ( i )->mean();
+	}
 	return mu;
 }
@@ -98,20 +109,22 @@
 	// non-central moment
 	vec mom2 = zeros ( dim );
-	for ( int i = 0;i < w.length();i++ ) {
+	for ( int i = 0; i < w.length(); i++ ) {
 		// pow is overloaded, we have to use another approach
-		mom2 += w ( i ) * (Coms(i)->variance() + elem_mult ( Coms(i)->mean(), Coms(i)->mean() )); 
+		mom2 += w ( i ) * ( Coms ( i )->variance() + elem_mult ( Coms ( i )->mean(), Coms ( i )->mean() ) );
 	}
 	// central moment
 	// pow is overloaded, we have to use another approach
-	return mom2 - elem_mult( mean(), mean() );
-}
-
-emix* egiwmix::marginal(const RV &rv) const{
-	it_assert_debug(isnamed(), "rvs are not assigned");
-			
-	Array<epdf*> Cn(Coms.length());
-	for(int i=0;i<Coms.length();i++){Cn(i)=Coms(i)->marginal(rv);}
+	return mom2 - elem_mult ( mean(), mean() );
+}
+
+emix* egiwmix::marginal ( const RV &rv ) const {
+	it_assert_debug ( isnamed(), "rvs are not assigned" );
+
+	Array<epdf*> Cn ( Coms.length() );
+	for ( int i = 0; i < Coms.length(); i++ ) {
+		Cn ( i ) = Coms ( i )->marginal ( rv );
+	}
 	emix* tmp = new emix();
-	tmp->set_parameters(w,Cn,false);
+	tmp->set_parameters ( w, Cn, false );
 	tmp->ownComs();
 	return tmp;
@@ -124,39 +137,40 @@
 
 	double sumVecCommon;			      	// common part for many terms in eq.
-	int len = w.length();				// no. of mix components	
-	int dimLS = Coms(1)->_V()._D().length() - 1; 	// dim of LS
-	vec vecNu(len);					// vector of dfms of components
-	vec vecD(len);					// vector of LS reminders of comps.
-	vec vecCommon(len);				// vector of common parts
+	int len = w.length();				// no. of mix components
+	int dimLS = Coms ( 1 )->_V()._D().length() - 1; 	// dim of LS
+	vec vecNu ( len );					// vector of dfms of components
+	vec vecD ( len );					// vector of LS reminders of comps.
+	vec vecCommon ( len );				// vector of common parts
 	mat matVecsTheta;				// matrix which rows are theta vects.
 
 	// fill in the vectors vecNu, vecD and matVecsTheta
-	for ( int i=0; i<len; i++ ) {
-		vecNu.shift_left( Coms(i)->_nu() );
-		vecD.shift_left( Coms(i)->_V()._D()(0) );
-		matVecsTheta.append_row( Coms(i)->est_theta()  );
+	for ( int i = 0; i < len; i++ ) {
+		vecNu.shift_left ( Coms ( i )->_nu() );
+		vecD.shift_left ( Coms ( i )->_V()._D() ( 0 ) );
+		matVecsTheta.append_row ( Coms ( i )->est_theta() );
 	}
 
 	// calculate the common parts and their sum
-	vecCommon = elem_mult ( w, elem_div(vecNu, vecD) );
-	sumVecCommon = sum(vecCommon);
+	vecCommon = elem_mult ( w, elem_div ( vecNu, vecD ) );
+	sumVecCommon = sum ( vecCommon );
 
 	// LS estimator of theta
-	vec aprEstTheta(dimLS);  aprEstTheta.zeros();
-	for ( int i=0; i<len; i++ ) {
-		aprEstTheta +=  matVecsTheta.get_row( i ) * vecCommon ( i );
+	vec aprEstTheta ( dimLS );
+	aprEstTheta.zeros();
+	for ( int i = 0; i < len; i++ ) {
+		aprEstTheta +=  matVecsTheta.get_row ( i ) * vecCommon ( i );
 	}
 	aprEstTheta /= sumVecCommon;
-	
-	
+
+
 	// LS estimator of dfm
 	double aprNu;
-	double A = log( sumVecCommon );		// Term 'A' in equation
-
-	for ( int i=0; i<len; i++ ) {
-		A += w(i) * ( log( vecD(i) ) - psi( 0.5 * vecNu(i) ) );
-	}
-
-	aprNu = ( 1 + sqrt(1 + 2 * (A - LOG2)/3 ) ) / ( 2 * (A - LOG2) );
+	double A = log ( sumVecCommon );		// Term 'A' in equation
+
+	for ( int i = 0; i < len; i++ ) {
+		A += w ( i ) * ( log ( vecD ( i ) ) - psi ( 0.5 * vecNu ( i ) ) );
+	}
+
+	aprNu = ( 1 + sqrt ( 1 + 2 * ( A - LOG2 ) / 3 ) ) / ( 2 * ( A - LOG2 ) );
 
 
@@ -167,20 +181,20 @@
 	// the following code is very numerically sensitive, thus
 	// we have to eliminate decompositions etc. as much as possible
-	mat aprC = zeros(dimLS, dimLS);
-	for ( int i=0; i<len; i++ ) {
-		aprC += Coms(i)->est_theta_cov().to_mat() * w(i); 
-		vec tmp = ( matVecsTheta.get_row(i) - aprEstTheta );
-		aprC += vecCommon(i) * outer_product( tmp, tmp);
+	mat aprC = zeros ( dimLS, dimLS );
+	for ( int i = 0; i < len; i++ ) {
+		aprC += Coms ( i )->est_theta_cov().to_mat() * w ( i );
+		vec tmp = ( matVecsTheta.get_row ( i ) - aprEstTheta );
+		aprC += vecCommon ( i ) * outer_product ( tmp, tmp );
 	}
 
 	// Construct GiW pdf :: BEGIN
-	ldmat aprCinv ( inv(aprC) );
-	vec D = concat( aprD, aprCinv._D() );
-	mat L = eye(dimLS+1);
-	L.set_submatrix(1,0, aprCinv._L() * aprEstTheta);
-	L.set_submatrix(1,1, aprCinv._L());
-	ldmat aprLD (L, D);
-
-	egiw* aprgiw = new egiw(1, aprLD, aprNu);
+	ldmat aprCinv ( inv ( aprC ) );
+	vec D = concat ( aprD, aprCinv._D() );
+	mat L = eye ( dimLS + 1 );
+	L.set_submatrix ( 1, 0, aprCinv._L() * aprEstTheta );
+	L.set_submatrix ( 1, 1, aprCinv._L() );
+	ldmat aprLD ( L, D );
+
+	egiw* aprgiw = new egiw ( 1, aprLD, aprNu );
 	return aprgiw;
 };
Index: library/bdm/stat/merger.h
===================================================================
--- library/bdm/stat/merger.h (revision 471)
+++ library/bdm/stat/merger.h (revision 477)
@@ -17,6 +17,5 @@
 #include "../estim/mixtures.h"
 
-namespace bdm
-{
+namespace bdm {
 using std::string;
 
@@ -43,330 +42,345 @@
 */
 
-class merger_base : public compositepdf, public epdf
-{
-	protected:
-		//! Data link for each mpdf in mpdfs
-		Array<datalink_m2e*> dls;
-		//! Array of rvs that are not modelled by mpdfs at all, \f$ z_i \f$
-		Array<RV> rvzs;
-		//! Data Links for extension \f$ f(z_i|x_i,y_i) \f$
-		Array<datalink_m2e*> zdls;
-		//! number of support points
-		int Npoints;
-		//! number of sources
-		int Nsources;
-
-		//! switch of the methoh used for merging
-		MERGER_METHOD METHOD;
-		//! Default for METHOD
-		static const MERGER_METHOD DFLT_METHOD;
-		
-		//!Prior on the log-normal merging model
-		double beta;
-		//! default for beta
-		static const double DFLT_beta;
-		
-		//! Projection to empirical density (could also be piece-wise linear)
-		eEmp eSmp;
-
-		//! debug or not debug
-		bool DBG;
-
-		//! debugging file
-		it_file* dbg_file;
-	public:
-		//! \name Constructors
-		//! @{
-
-		//!Empty constructor
-		merger_base () : compositepdf() {DBG = false;dbg_file = NULL;}
-
-		//!Constructor from sources
-		merger_base (const Array<mpdf*> &S, bool own=false);
-
-		//! Function setting the main internal structures
-		void set_sources (const Array<mpdf*> &Sources, bool own) {
-			compositepdf::set_elements (Sources,own);
-			Nsources=mpdfs.length();
-			//set sizes
-			dls.set_size (Sources.length());
-			rvzs.set_size (Sources.length());
-			zdls.set_size (Sources.length());
-
-			rv = getrv (/* checkoverlap = */ false);
-			RV rvc; setrvc (rv, rvc);  // Extend rv by rvc!
-			// join rv and rvc - see descriprion
-			rv.add (rvc);
-			// get dimension
-			dim = rv._dsize();
-
-			// create links between sources and common rv
-			RV xytmp;
-			for (int i = 0;i < mpdfs.length();i++) {
-				//Establich connection between mpdfs and merger
-				dls (i) = new datalink_m2e;
-				dls (i)->set_connection (mpdfs (i)->_rv(), mpdfs (i)->_rvc(), rv);
-
-				// find out what is missing in each mpdf
-				xytmp = mpdfs (i)->_rv();
-				xytmp.add (mpdfs (i)->_rvc());
-				// z_i = common_rv-xy
-				rvzs (i) = rv.subt (xytmp);
-				//establish connection between extension (z_i|x,y)s and common rv
-				zdls (i) = new datalink_m2e; zdls (i)->set_connection (rvzs (i), xytmp, rv) ;
-			};
-		}
-		//! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. Same number of points, \c dimsize, in each dimension. 
-		void set_support (const Array<vec> &XYZ, const int dimsize) {
-			set_support(XYZ,dimsize*ones_i(XYZ.length()));
-		}
-		//! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. \c gridsize specifies number of points is each dimension.
-		void set_support (const Array<vec> &XYZ, const ivec &gridsize) {
-			int dim = XYZ.length();  //check with internal dim!!
-			Npoints = prod (gridsize); 
-			eSmp.set_parameters (Npoints, false);
-			Array<vec> &samples = eSmp._samples();
-			eSmp._w() = ones (Npoints) / Npoints; //unifrom size of bins
-			//set samples
-			ivec ind = zeros_i (dim);      //indeces of dimensions in for cycle;
-			vec smpi (dim);            // ith sample
-			vec steps =zeros(dim);            // ith sample
-			// first corner
-			for (int j = 0; j < dim; j++) {
-				smpi (j) = XYZ (j) (0); /* beginning of the interval*/ 
-				it_assert(gridsize(j)!=0.0,"Zeros in gridsize!");
-				steps (j) = ( XYZ(j)(1)-smpi(j) )/gridsize(j);
-			}
-			// fill samples
-			for (int i = 0; i < Npoints; i++) {
-				// copy 
-				samples(i) = smpi; 
-				// go through all dimensions
-				for (int j = 0;j < dim;j++) {
-					if (ind (j) == gridsize (j) - 1) { //j-th index is full
-						ind (j) = 0; //shift back
-						smpi(j) = XYZ(j)(0);
-						
-						if (i<Npoints-1) {
-							ind (j + 1) ++; //increase the next dimension;
-							smpi(j+1) += steps(j+1);
-							break;
-						}
-						
-					} else {
-						ind (j) ++; 
-						smpi(j) +=steps(j);
+class merger_base : public compositepdf, public epdf {
+protected:
+	//! Data link for each mpdf in mpdfs
+	Array<datalink_m2e*> dls;
+	//! Array of rvs that are not modelled by mpdfs at all, \f$ z_i \f$
+	Array<RV> rvzs;
+	//! Data Links for extension \f$ f(z_i|x_i,y_i) \f$
+	Array<datalink_m2e*> zdls;
+	//! number of support points
+	int Npoints;
+	//! number of sources
+	int Nsources;
+
+	//! switch of the methoh used for merging
+	MERGER_METHOD METHOD;
+	//! Default for METHOD
+	static const MERGER_METHOD DFLT_METHOD;
+
+	//!Prior on the log-normal merging model
+	double beta;
+	//! default for beta
+	static const double DFLT_beta;
+
+	//! Projection to empirical density (could also be piece-wise linear)
+	eEmp eSmp;
+
+	//! debug or not debug
+	bool DBG;
+
+	//! debugging file
+	it_file* dbg_file;
+public:
+	//! \name Constructors
+	//! @{
+
+	//!Empty constructor
+	merger_base () : compositepdf() {
+		DBG = false;
+		dbg_file = NULL;
+	}
+
+	//!Constructor from sources
+	merger_base ( const Array<mpdf*> &S, bool own = false );
+
+	//! Function setting the main internal structures
+	void set_sources ( const Array<mpdf*> &Sources, bool own ) {
+		compositepdf::set_elements ( Sources, own );
+		Nsources = mpdfs.length();
+		//set sizes
+		dls.set_size ( Sources.length() );
+		rvzs.set_size ( Sources.length() );
+		zdls.set_size ( Sources.length() );
+
+		rv = getrv ( /* checkoverlap = */ false );
+		RV rvc;
+		setrvc ( rv, rvc );  // Extend rv by rvc!
+		// join rv and rvc - see descriprion
+		rv.add ( rvc );
+		// get dimension
+		dim = rv._dsize();
+
+		// create links between sources and common rv
+		RV xytmp;
+		for ( int i = 0; i < mpdfs.length(); i++ ) {
+			//Establich connection between mpdfs and merger
+			dls ( i ) = new datalink_m2e;
+			dls ( i )->set_connection ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv );
+
+			// find out what is missing in each mpdf
+			xytmp = mpdfs ( i )->_rv();
+			xytmp.add ( mpdfs ( i )->_rvc() );
+			// z_i = common_rv-xy
+			rvzs ( i ) = rv.subt ( xytmp );
+			//establish connection between extension (z_i|x,y)s and common rv
+			zdls ( i ) = new datalink_m2e;
+			zdls ( i )->set_connection ( rvzs ( i ), xytmp, rv ) ;
+		};
+	}
+	//! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. Same number of points, \c dimsize, in each dimension.
+	void set_support ( const Array<vec> &XYZ, const int dimsize ) {
+		set_support ( XYZ, dimsize*ones_i ( XYZ.length() ) );
+	}
+	//! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. \c gridsize specifies number of points is each dimension.
+	void set_support ( const Array<vec> &XYZ, const ivec &gridsize ) {
+		int dim = XYZ.length();  //check with internal dim!!
+		Npoints = prod ( gridsize );
+		eSmp.set_parameters ( Npoints, false );
+		Array<vec> &samples = eSmp._samples();
+		eSmp._w() = ones ( Npoints ) / Npoints; //unifrom size of bins
+		//set samples
+		ivec ind = zeros_i ( dim );    //indeces of dimensions in for cycle;
+		vec smpi ( dim );          // ith sample
+		vec steps = zeros ( dim );        // ith sample
+		// first corner
+		for ( int j = 0; j < dim; j++ ) {
+			smpi ( j ) = XYZ ( j ) ( 0 ); /* beginning of the interval*/
+			it_assert ( gridsize ( j ) != 0.0, "Zeros in gridsize!" );
+			steps ( j ) = ( XYZ ( j ) ( 1 ) - smpi ( j ) ) / gridsize ( j );
+		}
+		// fill samples
+		for ( int i = 0; i < Npoints; i++ ) {
+			// copy
+			samples ( i ) = smpi;
+			// go through all dimensions
+			for ( int j = 0; j < dim; j++ ) {
+				if ( ind ( j ) == gridsize ( j ) - 1 ) { //j-th index is full
+					ind ( j ) = 0; //shift back
+					smpi ( j ) = XYZ ( j ) ( 0 );
+
+					if ( i < Npoints - 1 ) {
+						ind ( j + 1 ) ++; //increase the next dimension;
+						smpi ( j + 1 ) += steps ( j + 1 );
 						break;
 					}
+
+				} else {
+					ind ( j ) ++;
+					smpi ( j ) += steps ( j );
+					break;
 				}
 			}
 		}
-		//! set debug file
-		void set_debug_file (const string fname) {
-			if (DBG) delete dbg_file;
-			dbg_file = new it_file (fname);
-			if (dbg_file) DBG = true;
-		}
-		//! Set internal parameters used in approximation
-		void set_method (MERGER_METHOD MTH=DFLT_METHOD, double beta0 = DFLT_beta) {
-			METHOD = MTH;
-			beta = beta0;
-		}
-		//! Set support points from a pdf by drawing N samples
-		void set_support (const epdf &overall, int N) {
-			eSmp.set_statistics (&overall, N);
-			Npoints = N;
-		}
-
-		//! Destructor
-		virtual ~merger_base() {
-			for (int i = 0; i < Nsources; i++) {
-				delete dls (i);
-				delete zdls (i);
-			}
-			if (DBG) delete dbg_file;
-		};
-		//!@}
-
-		//! \name Mathematical operations
-		//!@{
-
-		//!Merge given sources in given points
-		virtual void merge () {
-			validate();
-
-			//check if sources overlap:
-			bool OK = true;
-			for (int i = 0;i < mpdfs.length(); i++) {
-				OK &= (rvzs (i)._dsize() == 0); // z_i is empty
-				OK &= (mpdfs (i)->_rvc()._dsize() == 0); // y_i is empty
-			}
-
-			if (OK) {
-				mat lW = zeros (mpdfs.length(), eSmp._w().length());
-
-				vec emptyvec (0);
-				for (int i = 0; i < mpdfs.length(); i++) {
-					for (int j = 0; j < eSmp._w().length(); j++) {
-						lW (i, j) = mpdfs (i)->evallogcond (eSmp._samples() (j), emptyvec);
-					}
-				}
-
-				vec w_nn=merge_points (lW);
-				vec wtmp = exp (w_nn-max(w_nn));
-				//renormalize
-				eSmp._w() = wtmp / sum (wtmp);
-			} else {
-				it_error ("Sources are not compatible - use merger_mix");
-			}
-		};
-
-
-		//! Merge log-likelihood values in points using method specified by parameter METHOD
-		vec merge_points (mat &lW);
-
-
-		//! sample from merged density
-//! weight w is a
-		vec mean() const {
-			const Vec<double> &w = eSmp._w();
-			const Array<vec> &S = eSmp._samples();
-			vec tmp = zeros (dim);
-			for (int i = 0; i < Npoints; i++) {
-				tmp += w (i) * S (i);
-			}
-			return tmp;
-		}
-		mat covariance() const {
-			const vec &w = eSmp._w();
-			const Array<vec> &S = eSmp._samples();
-
-			vec mea = mean();
-
-// 			cout << sum (w) << "," << w*w << endl;
-
-			mat Tmp = zeros (dim, dim);
-			for (int i = 0; i < Npoints; i++) {
-				Tmp += w (i) * outer_product (S (i), S (i));
-			}
-			return Tmp -outer_product (mea, mea);
-		}
-		vec variance() const {
-			const vec &w = eSmp._w();
-			const Array<vec> &S = eSmp._samples();
-
-			vec tmp = zeros (dim);
-			for (int i = 0; i < Nsources; i++) {
-				tmp += w (i) * pow (S (i), 2);
-			}
-			return tmp -pow (mean(), 2);
-		}
-		//!@}
-
-		//! \name Access to attributes
-		//! @{
-
-		//! Access function
-		eEmp& _Smp() {return eSmp;}
-
-		//! load from setting
-		void from_setting (const Setting& set) {
-			// get support
-			// find which method to use
-			string meth_str;
-			UI::get<string> (meth_str, set, "method", UI::compulsory);
-			if (!strcmp (meth_str.c_str(), "arithmetic"))
-				set_method (ARITHMETIC);
-			else {
-				if (!strcmp (meth_str.c_str(), "geometric"))
-					set_method (GEOMETRIC);
-				else if (!strcmp (meth_str.c_str(), "lognormal")) {
-					set_method (LOGNORMAL);
-					set.lookupValue( "beta",beta);
+	}
+	//! set debug file
+	void set_debug_file ( const string fname ) {
+		if ( DBG ) delete dbg_file;
+		dbg_file = new it_file ( fname );
+		if ( dbg_file ) DBG = true;
+	}
+	//! Set internal parameters used in approximation
+	void set_method ( MERGER_METHOD MTH = DFLT_METHOD, double beta0 = DFLT_beta ) {
+		METHOD = MTH;
+		beta = beta0;
+	}
+	//! Set support points from a pdf by drawing N samples
+	void set_support ( const epdf &overall, int N ) {
+		eSmp.set_statistics ( &overall, N );
+		Npoints = N;
+	}
+
+	//! Destructor
+	virtual ~merger_base() {
+		for ( int i = 0; i < Nsources; i++ ) {
+			delete dls ( i );
+			delete zdls ( i );
+		}
+		if ( DBG ) delete dbg_file;
+	};
+	//!@}
+
+	//! \name Mathematical operations
+	//!@{
+
+	//!Merge given sources in given points
+	virtual void merge () {
+		validate();
+
+		//check if sources overlap:
+		bool OK = true;
+		for ( int i = 0; i < mpdfs.length(); i++ ) {
+			OK &= ( rvzs ( i )._dsize() == 0 ); // z_i is empty
+			OK &= ( mpdfs ( i )->_rvc()._dsize() == 0 ); // y_i is empty
+		}
+
+		if ( OK ) {
+			mat lW = zeros ( mpdfs.length(), eSmp._w().length() );
+
+			vec emptyvec ( 0 );
+			for ( int i = 0; i < mpdfs.length(); i++ ) {
+				for ( int j = 0; j < eSmp._w().length(); j++ ) {
+					lW ( i, j ) = mpdfs ( i )->evallogcond ( eSmp._samples() ( j ), emptyvec );
 				}
 			}
-			string dbg_file;
-			if (UI::get(dbg_file, set, "dbg_file"))
-				set_debug_file(dbg_file);
-			//validate() - not used
-		}
-
-		void validate() {
-			it_assert (eSmp._w().length() > 0, "Empty support, use set_support().");
-			it_assert (dim == eSmp._samples() (0).length(), "Support points and rv are not compatible!");
-			it_assert (isnamed(),"mergers must be named");
-		}
-		//!@}
+
+			vec w_nn = merge_points ( lW );
+			vec wtmp = exp ( w_nn - max ( w_nn ) );
+			//renormalize
+			eSmp._w() = wtmp / sum ( wtmp );
+		} else {
+			it_error ( "Sources are not compatible - use merger_mix" );
+		}
+	};
+
+
+	//! Merge log-likelihood values in points using method specified by parameter METHOD
+	vec merge_points ( mat &lW );
+
+
+	//! sample from merged density
+//! weight w is a
+	vec mean() const {
+		const Vec<double> &w = eSmp._w();
+		const Array<vec> &S = eSmp._samples();
+		vec tmp = zeros ( dim );
+		for ( int i = 0; i < Npoints; i++ ) {
+			tmp += w ( i ) * S ( i );
+		}
+		return tmp;
+	}
+	mat covariance() const {
+		const vec &w = eSmp._w();
+		const Array<vec> &S = eSmp._samples();
+
+		vec mea = mean();
+
+// 			cout << sum (w) << "," << w*w << endl;
+
+		mat Tmp = zeros ( dim, dim );
+		for ( int i = 0; i < Npoints; i++ ) {
+			Tmp += w ( i ) * outer_product ( S ( i ), S ( i ) );
+		}
+		return Tmp - outer_product ( mea, mea );
+	}
+	vec variance() const {
+		const vec &w = eSmp._w();
+		const Array<vec> &S = eSmp._samples();
+
+		vec tmp = zeros ( dim );
+		for ( int i = 0; i < Nsources; i++ ) {
+			tmp += w ( i ) * pow ( S ( i ), 2 );
+		}
+		return tmp - pow ( mean(), 2 );
+	}
+	//!@}
+
+	//! \name Access to attributes
+	//! @{
+
+	//! Access function
+	eEmp& _Smp() {
+		return eSmp;
+	}
+
+	//! load from setting
+	void from_setting ( const Setting& set ) {
+		// get support
+		// find which method to use
+		string meth_str;
+		UI::get<string> ( meth_str, set, "method", UI::compulsory );
+		if ( !strcmp ( meth_str.c_str(), "arithmetic" ) )
+			set_method ( ARITHMETIC );
+		else {
+			if ( !strcmp ( meth_str.c_str(), "geometric" ) )
+				set_method ( GEOMETRIC );
+			else if ( !strcmp ( meth_str.c_str(), "lognormal" ) ) {
+				set_method ( LOGNORMAL );
+				set.lookupValue ( "beta", beta );
+			}
+		}
+		string dbg_file;
+		if ( UI::get ( dbg_file, set, "dbg_file" ) )
+			set_debug_file ( dbg_file );
+		//validate() - not used
+	}
+
+	void validate() {
+		it_assert ( eSmp._w().length() > 0, "Empty support, use set_support()." );
+		it_assert ( dim == eSmp._samples() ( 0 ).length(), "Support points and rv are not compatible!" );
+		it_assert ( isnamed(), "mergers must be named" );
+	}
+	//!@}
 };
-UIREGISTER(merger_base);
-
-class merger_mix : public merger_base
-{
-	protected:
-		//!Internal mixture of EF models
-		MixEF Mix;
-		//!Number of components in a mixture
-		int Ncoms;
-		//! coefficient of resampling [0,1]
-		double effss_coef;
-		//! stop after niter iterations
-		int stop_niter;
-		
-		//! default value for Ncoms
-		static const int DFLT_Ncoms;
-		//! default value for efss_coef;
-		static const double DFLT_effss_coef;
-
-	public:
-		//!\name Constructors
-		//!@{
-		merger_mix () {};
-		merger_mix (const Array<mpdf*> &S,bool own=false) {set_sources (S,own);};
-		//! Set sources and prepare all internal structures
-		void set_sources (const Array<mpdf*> &S, bool own) {
-			merger_base::set_sources (S,own);
-			Nsources = S.length();
-		}
-		//! Set internal parameters used in approximation
-		void set_parameters (int Ncoms0 = DFLT_Ncoms, double effss_coef0 = DFLT_effss_coef) {
-			Ncoms = Ncoms0;
-			effss_coef = effss_coef0;
-		}
-		//!@}
-
-		//! \name Mathematical operations
-		//!@{
-
-		//!Merge values using mixture approximation
-		void merge ();
-
-		//! sample from the approximating mixture
-		vec sample () const { return Mix.posterior().sample();}
-		//! loglikelihood computed on mixture models
-		double evallog (const vec &dt) const {
-			vec dtf = ones (dt.length() + 1);
-			dtf.set_subvector (0, dt);
-			return Mix.logpred (dtf);
-		}
-		//!@}
-
-		//!\name Access functions
-		//!@{
+UIREGISTER ( merger_base );
+
+class merger_mix : public merger_base {
+protected:
+	//!Internal mixture of EF models
+	MixEF Mix;
+	//!Number of components in a mixture
+	int Ncoms;
+	//! coefficient of resampling [0,1]
+	double effss_coef;
+	//! stop after niter iterations
+	int stop_niter;
+
+	//! default value for Ncoms
+	static const int DFLT_Ncoms;
+	//! default value for efss_coef;
+	static const double DFLT_effss_coef;
+
+public:
+	//!\name Constructors
+	//!@{
+	merger_mix () {};
+	merger_mix ( const Array<mpdf*> &S, bool own = false ) {
+		set_sources ( S, own );
+	};
+	//! Set sources and prepare all internal structures
+	void set_sources ( const Array<mpdf*> &S, bool own ) {
+		merger_base::set_sources ( S, own );
+		Nsources = S.length();
+	}
+	//! Set internal parameters used in approximation
+	void set_parameters ( int Ncoms0 = DFLT_Ncoms, double effss_coef0 = DFLT_effss_coef ) {
+		Ncoms = Ncoms0;
+		effss_coef = effss_coef0;
+	}
+	//!@}
+
+	//! \name Mathematical operations
+	//!@{
+
+	//!Merge values using mixture approximation
+	void merge ();
+
+	//! sample from the approximating mixture
+	vec sample () const {
+		return Mix.posterior().sample();
+	}
+	//! loglikelihood computed on mixture models
+	double evallog ( const vec &dt ) const {
+		vec dtf = ones ( dt.length() + 1 );
+		dtf.set_subvector ( 0, dt );
+		return Mix.logpred ( dtf );
+	}
+	//!@}
+
+	//!\name Access functions
+	//!@{
 //! Access function
-		MixEF& _Mix() {return Mix;}
-		//! Access function
-		emix* proposal() {emix* tmp = Mix.epredictor(); tmp->set_rv (rv); return tmp;}
-		//! from_settings
-		void from_setting(const Setting& set){
-			merger_base::from_setting(set);
-			set.lookupValue("ncoms",Ncoms);
-			set.lookupValue("effss_coef",effss_coef);
-			set.lookupValue("stop_niter",stop_niter);
-		}
-		
-		//! @}
+	MixEF& _Mix() {
+		return Mix;
+	}
+	//! Access function
+	emix* proposal() {
+		emix* tmp = Mix.epredictor();
+		tmp->set_rv ( rv );
+		return tmp;
+	}
+	//! from_settings
+	void from_setting ( const Setting& set ) {
+		merger_base::from_setting ( set );
+		set.lookupValue ( "ncoms", Ncoms );
+		set.lookupValue ( "effss_coef", effss_coef );
+		set.lookupValue ( "stop_niter", stop_niter );
+	}
+
+	//! @}
 
 };
-UIREGISTER(merger_mix);
+UIREGISTER ( merger_mix );
 
 }
Index: library/bdm/stat/emix.h
===================================================================
--- library/bdm/stat/emix.h (revision 471)
+++ library/bdm/stat/emix.h (revision 477)
@@ -14,5 +14,5 @@
 #define EMIX_H
 
-#define LOG2  0.69314718055995	
+#define LOG2  0.69314718055995
 
 #include "../shared_ptr.h"
@@ -48,33 +48,45 @@
 	//!Default constructor. By default, the given epdf is not copied!
 	//! It is assumed that this function will be used only temporarily.
-	mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( ), dl ( ) {
+	mratio ( const epdf* nom0, const RV &rv, bool copy = false ) : mpdf ( ), dl ( ) {
 		// adjust rv and rvc
 		rvc = nom0->_rv().subt ( rv );
 		dimc = rvc._dsize();
-		set_ep(shared_ptr<epdf>(new epdf));
-		e()->set_parameters(rv._dsize());
-		e()->set_rv(rv);
-		
+		set_ep ( shared_ptr<epdf> ( new epdf ) );
+		e()->set_parameters ( rv._dsize() );
+		e()->set_rv ( rv );
+
 		//prepare data structures
-		if ( copy ) {it_error ( "todo" ); destroynom=true; }
-		else { nom = nom0; destroynom = false; }
-		it_assert_debug ( rvc.length() >0,"Makes no sense to use this object!" );		
-		
+		if ( copy ) {
+			it_error ( "todo" );
+			destroynom = true;
+		} else {
+			nom = nom0;
+			destroynom = false;
+		}
+		it_assert_debug ( rvc.length() > 0, "Makes no sense to use this object!" );
+
 		// build denominator
 		den = nom->marginal ( rvc );
-		dl.set_connection(rv,rvc,nom0->_rv());
+		dl.set_connection ( rv, rvc, nom0->_rv() );
 	};
 	double evallogcond ( const vec &val, const vec &cond ) {
 		double tmp;
 		vec nom_val ( e()->dimension() + dimc );
-		dl.pushup_cond ( nom_val,val,cond );
+		dl.pushup_cond ( nom_val, val, cond );
 		tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
-		it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+		it_assert_debug ( std::isfinite ( tmp ), "Infinite value" );
 		return tmp;
 	}
 	//! Object takes ownership of nom and will destroy it
-	void ownnom() {destroynom=true;}
+	void ownnom() {
+		destroynom = true;
+	}
 	//! Default destructor
-	~mratio() {delete den; if ( destroynom ) {delete nom;}}
+	~mratio() {
+		delete den;
+		if ( destroynom ) {
+			delete nom;
+		}
+	}
 };
 
@@ -102,10 +114,13 @@
 	//! Set weights \c w and components \c Coms
 	//!By default Coms are copied inside. Parameter \c copy can be set to false if Coms live externally. Use method ownComs() if Coms should be destroyed by the destructor.
-	void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy=false );
+	void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy = false );
 
 	vec sample() const;
 	vec mean() const {
-		int i; vec mu = zeros ( dim );
-		for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
+		int i;
+		vec mu = zeros ( dim );
+		for ( i = 0; i < w.length(); i++ ) {
+			mu += w ( i ) * Coms ( i )->mean();
+		}
 		return mu;
 	}
@@ -113,21 +128,27 @@
 		//non-central moment
 		vec mom2 = zeros ( dim );
-		for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * (Coms(i)->variance() + pow ( Coms ( i )->mean(),2 )); }
+		for ( int i = 0; i < w.length(); i++ ) {
+			mom2 += w ( i ) * ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) );
+		}
 		//central moment
-		return mom2-pow ( mean(),2 );
+		return mom2 - pow ( mean(), 2 );
 	}
 	double evallog ( const vec &val ) const {
 		int i;
 		double sum = 0.0;
-		for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );}
-		if ( sum==0.0 ) {sum=std::numeric_limits<double>::epsilon();}
-		double tmp=log ( sum );
-		it_assert_debug ( std::isfinite ( tmp ),"Infinite" );
+		for ( i = 0; i < w.length(); i++ ) {
+			sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );
+		}
+		if ( sum == 0.0 ) {
+			sum = std::numeric_limits<double>::epsilon();
+		}
+		double tmp = log ( sum );
+		it_assert_debug ( std::isfinite ( tmp ), "Infinite" );
 		return tmp;
 	};
 	vec evallog_m ( const mat &Val ) const {
-		vec x=zeros ( Val.cols() );
+		vec x = zeros ( Val.cols() );
 		for ( int i = 0; i < w.length(); i++ ) {
-			x+= w ( i ) *exp ( Coms ( i )->evallog_m ( Val ) );
+			x += w ( i ) * exp ( Coms ( i )->evallog_m ( Val ) );
 		}
 		return log ( x );
@@ -147,14 +168,28 @@
 //Access methods
 	//! returns a pointer to the internal mean value. Use with Care!
-	vec& _w() {return w;}
-	virtual ~emix() {if ( destroyComs ) {for ( int i=0;i<Coms.length();i++ ) {delete Coms ( i );}}}
+	vec& _w() {
+		return w;
+	}
+	virtual ~emix() {
+		if ( destroyComs ) {
+			for ( int i = 0; i < Coms.length(); i++ ) {
+				delete Coms ( i );
+			}
+		}
+	}
 	//! Auxiliary function for taking ownership of the Coms()
-	void ownComs() {destroyComs=true;}
+	void ownComs() {
+		destroyComs = true;
+	}
 
 	//!access function
-	epdf* _Coms ( int i ) {return Coms ( i );}
-	void set_rv(const RV &rv){
-		epdf::set_rv(rv);
-		for(int i=0;i<Coms.length();i++){Coms(i)->set_rv(rv);}
+	epdf* _Coms ( int i ) {
+		return Coms ( i );
+	}
+	void set_rv ( const RV &rv ) {
+		epdf::set_rv ( rv );
+		for ( int i = 0; i < Coms.length(); i++ ) {
+			Coms ( i )->set_rv ( rv );
+		}
 	}
 };
@@ -179,5 +214,5 @@
 	//! Set weights \c w and components \c Coms
 	//!By default Coms are copied inside. Parameter \c copy can be set to false if Coms live externally. Use method ownComs() if Coms should be destroyed by the destructor.
-	void set_parameters ( const vec &w, const Array<egiw*> &Coms, bool copy=false );
+	void set_parameters ( const vec &w, const Array<egiw*> &Coms, bool copy = false );
 
 	//!return expected value
@@ -188,25 +223,43 @@
 
 	//!return the expected variance
- 	vec variance() const;
+	vec variance() const;
 
 	// TODO!!! Defined to follow ANSI and/or for future development
 	void mean_mat ( mat &M, mat&R ) const {};
-	double evallog_nn ( const vec &val ) const {return 0;};
-	double lognc () const {return 0;};
+	double evallog_nn ( const vec &val ) const {
+		return 0;
+	};
+	double lognc () const {
+		return 0;
+	};
 	emix* marginal ( const RV &rv ) const;
 
 //Access methods
 	//! returns a pointer to the internal mean value. Use with Care!
-	vec& _w() {return w;}
-	virtual ~egiwmix() {if ( destroyComs ) {for ( int i=0;i<Coms.length();i++ ) {delete Coms ( i );}}}
+	vec& _w() {
+		return w;
+	}
+	virtual ~egiwmix() {
+		if ( destroyComs ) {
+			for ( int i = 0; i < Coms.length(); i++ ) {
+				delete Coms ( i );
+			}
+		}
+	}
 	//! Auxiliary function for taking ownership of the Coms()
-	void ownComs() {destroyComs=true;}
+	void ownComs() {
+		destroyComs = true;
+	}
 
 	//!access function
-	egiw* _Coms ( int i ) {return Coms ( i );}
-
-	void set_rv(const RV &rv){
-		egiw::set_rv(rv);
-		for(int i=0;i<Coms.length();i++){Coms(i)->set_rv(rv);}
+	egiw* _Coms ( int i ) {
+		return Coms ( i );
+	}
+
+	void set_rv ( const RV &rv ) {
+		egiw::set_rv ( rv );
+		for ( int i = 0; i < Coms.length(); i++ ) {
+			Coms ( i )->set_rv ( rv );
+		}
 	}
 
@@ -236,29 +289,29 @@
 	/*!\brief Constructor from list of mFacs,
 	*/
-	mprod():dummy(new epdf()) { }
-	mprod (Array<mpdf*> mFacs):
-	    dummy(new epdf()) {
-	    set_elements(mFacs);
-	}
-
-	void set_elements(Array<mpdf*> mFacs , bool own=false) {
-		
-		compositepdf::set_elements(mFacs,own);
-		dls.set_size(mFacs.length());
-		epdfs.set_size(mFacs.length());
-				
-		set_ep(dummy);
-		RV rv=getrv ( true );
+	mprod() : dummy ( new epdf() ) { }
+	mprod ( Array<mpdf*> mFacs ) :
+			dummy ( new epdf() ) {
+		set_elements ( mFacs );
+	}
+
+	void set_elements ( Array<mpdf*> mFacs , bool own = false ) {
+
+		compositepdf::set_elements ( mFacs, own );
+		dls.set_size ( mFacs.length() );
+		epdfs.set_size ( mFacs.length() );
+
+		set_ep ( dummy );
+		RV rv = getrv ( true );
 		set_rv ( rv );
-		dummy->set_parameters(rv._dsize());
-		setrvc(e()->_rv(), rvc);
+		dummy->set_parameters ( rv._dsize() );
+		setrvc ( e()->_rv(), rvc );
 		// rv and rvc established = > we can link them with mpdfs
-		for ( int i = 0;i < mpdfs.length(); i++ ) {
+		for ( int i = 0; i < mpdfs.length(); i++ ) {
 			dls ( i ) = new datalink_m2m;
-			dls(i)->set_connection( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), _rv(), _rvc() );
-		}
-
-		for ( int i=0; i<mpdfs.length(); i++ ) {
-			epdfs(i) = mpdfs(i)->e();
+			dls ( i )->set_connection ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), _rv(), _rvc() );
+		}
+
+		for ( int i = 0; i < mpdfs.length(); i++ ) {
+			epdfs ( i ) = mpdfs ( i )->e();
 		}
 	};
@@ -267,5 +320,5 @@
 		int i;
 		double res = 0.0;
-		for ( i = mpdfs.length() - 1;i >= 0;i-- ) {
+		for ( i = mpdfs.length() - 1; i >= 0; i-- ) {
 			/*			if ( mpdfs(i)->_rvc().count() >0) {
 							mpdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
@@ -280,17 +333,17 @@
 		return res;
 	}
-	vec evallogcond_m(const mat &Dt, const vec &cond) {
-		vec tmp(Dt.cols());
-		for(int i=0;i<Dt.cols(); i++){
-			tmp(i) = evallogcond(Dt.get_col(i),cond);
-		}
-		return tmp;
-	};
-	vec evallogcond_m(const Array<vec> &Dt, const vec &cond) {
-		vec tmp(Dt.length());
-		for(int i=0;i<Dt.length(); i++){
-			tmp(i) = evallogcond(Dt(i),cond);
-		}
-		return tmp;		
+	vec evallogcond_m ( const mat &Dt, const vec &cond ) {
+		vec tmp ( Dt.cols() );
+		for ( int i = 0; i < Dt.cols(); i++ ) {
+			tmp ( i ) = evallogcond ( Dt.get_col ( i ), cond );
+		}
+		return tmp;
+	};
+	vec evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
+		vec tmp ( Dt.length() );
+		for ( int i = 0; i < Dt.length(); i++ ) {
+			tmp ( i ) = evallogcond ( Dt ( i ), cond );
+		}
+		return tmp;
 	};
 
@@ -298,11 +351,11 @@
 	//TODO smarter...
 	vec samplecond ( const vec &cond ) {
-	    //! Ugly hack to help to discover if mpfs are not in proper order. Correct solution = check that explicitely.
-	    vec smp= std::numeric_limits<double>::infinity() * ones ( e()->dimension() );
+		//! Ugly hack to help to discover if mpfs are not in proper order. Correct solution = check that explicitely.
+		vec smp = std::numeric_limits<double>::infinity() * ones ( e()->dimension() );
 		vec smpi;
 		// Hard assumption here!!! We are going backwards, to assure that samples that are needed from smp are already generated!
-		for ( int i = ( mpdfs.length() - 1 );i >= 0;i-- ) {
+		for ( int i = ( mpdfs.length() - 1 ); i >= 0; i-- ) {
 			if ( mpdfs ( i )->dimensionc() ) {
-				mpdfs ( i )->condition ( dls ( i )->get_cond ( smp ,cond ) ); // smp is val here!!
+				mpdfs ( i )->condition ( dls ( i )->get_cond ( smp , cond ) ); // smp is val here!!
 			}
 			smpi = epdfs ( i )->sample();
@@ -314,6 +367,8 @@
 	}
 	mat samplecond ( const vec &cond,  int N ) {
-		mat Smp ( dimension(),N );
-		for ( int i=0;i<N;i++ ) {Smp.set_col ( i,samplecond ( cond ) );}
+		mat Smp ( dimension(), N );
+		for ( int i = 0; i < N; i++ ) {
+			Smp.set_col ( i, samplecond ( cond ) );
+		}
 		return Smp;
 	}
@@ -327,12 +382,12 @@
 	//! \endcode
 	//!@}
-	void from_setting(const Setting &set){
+	void from_setting ( const Setting &set ) {
 		Array<mpdf*> Atmp; //temporary Array
-		UI::get(Atmp,set, "mpdfs", UI::compulsory);
-		set_elements(Atmp,true);
-	}
-	
+		UI::get ( Atmp, set, "mpdfs", UI::compulsory );
+		set_elements ( Atmp, true );
+	}
+
 };
-UIREGISTER(mprod);
+UIREGISTER ( mprod );
 
 //! Product of independent epdfs. For dependent pdfs, use mprod.
@@ -344,33 +399,34 @@
 	Array<datalink*> dls;
 public:
-	eprod () : epdfs ( 0 ),dls ( 0 ) {};
-	void set_parameters ( const Array<const epdf*> &epdfs0, bool named=true ) {
-		epdfs=epdfs0;//.set_length ( epdfs0.length() );
+	eprod () : epdfs ( 0 ), dls ( 0 ) {};
+	void set_parameters ( const Array<const epdf*> &epdfs0, bool named = true ) {
+		epdfs = epdfs0;//.set_length ( epdfs0.length() );
 		dls.set_length ( epdfs.length() );
 
-		bool independent=true;
+		bool independent = true;
 		if ( named ) {
-			for ( int i=0;i<epdfs.length();i++ ) {
-				independent=rv.add ( epdfs ( i )->_rv() );
-				it_assert_debug ( independent==true, "eprod:: given components are not independent." );
+			for ( int i = 0; i < epdfs.length(); i++ ) {
+				independent = rv.add ( epdfs ( i )->_rv() );
+				it_assert_debug ( independent == true, "eprod:: given components are not independent." );
 			}
-			dim=rv._dsize();
-		}
-		else {
-			dim =0; for ( int i=0;i<epdfs.length();i++ ) {
-				dim+=epdfs ( i )->dimension();
+			dim = rv._dsize();
+		} else {
+			dim = 0;
+			for ( int i = 0; i < epdfs.length(); i++ ) {
+				dim += epdfs ( i )->dimension();
 			}
 		}
 		//
-		int cumdim=0;
-		int dimi=0;
+		int cumdim = 0;
+		int dimi = 0;
 		int i;
-		for ( i=0;i<epdfs.length();i++ ) {
+		for ( i = 0; i < epdfs.length(); i++ ) {
 			dls ( i ) = new datalink;
-			if ( named ) {dls ( i )->set_connection ( epdfs ( i )->_rv() , rv );}
-			else {
+			if ( named ) {
+				dls ( i )->set_connection ( epdfs ( i )->_rv() , rv );
+			} else {
 				dimi = epdfs ( i )->dimension();
-				dls ( i )->set_connection ( dimi, dim, linspace ( cumdim,cumdim+dimi-1 ) );
-				cumdim+=dimi;
+				dls ( i )->set_connection ( dimi, dim, linspace ( cumdim, cumdim + dimi - 1 ) );
+				cumdim += dimi;
 			}
 		}
@@ -379,5 +435,5 @@
 	vec mean() const {
 		vec tmp ( dim );
-		for ( int i=0;i<epdfs.length();i++ ) {
+		for ( int i = 0; i < epdfs.length(); i++ ) {
 			vec pom = epdfs ( i )->mean();
 			dls ( i )->pushup ( tmp, pom );
@@ -387,13 +443,13 @@
 	vec variance() const {
 		vec tmp ( dim ); //second moment
-		for ( int i=0;i<epdfs.length();i++ ) {
+		for ( int i = 0; i < epdfs.length(); i++ ) {
 			vec pom = epdfs ( i )->mean();
-			dls ( i )->pushup ( tmp, pow ( pom,2 ) );
-		}
-		return tmp-pow ( mean(),2 );
+			dls ( i )->pushup ( tmp, pow ( pom, 2 ) );
+		}
+		return tmp - pow ( mean(), 2 );
 	}
 	vec sample() const {
 		vec tmp ( dim );
-		for ( int i=0;i<epdfs.length();i++ ) {
+		for ( int i = 0; i < epdfs.length(); i++ ) {
 			vec pom = epdfs ( i )->sample();
 			dls ( i )->pushup ( tmp, pom );
@@ -402,16 +458,23 @@
 	}
 	double evallog ( const vec &val ) const {
-		double tmp=0;
-		for ( int i=0;i<epdfs.length();i++ ) {
-			tmp+=epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
-		}
-		it_assert_debug ( std::isfinite ( tmp ),"Infinite" );
+		double tmp = 0;
+		for ( int i = 0; i < epdfs.length(); i++ ) {
+			tmp += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
+		}
+		it_assert_debug ( std::isfinite ( tmp ), "Infinite" );
 		return tmp;
 	}
 	//!access function
-	const epdf* operator () ( int i ) const {it_assert_debug ( i<epdfs.length(),"wrong index" );return epdfs ( i );}
+	const epdf* operator () ( int i ) const {
+		it_assert_debug ( i < epdfs.length(), "wrong index" );
+		return epdfs ( i );
+	}
 
 	//!Destructor
-	~eprod() {for ( int i=0;i<epdfs.length();i++ ) {delete dls ( i );}}
+	~eprod() {
+		for ( int i = 0; i < epdfs.length(); i++ ) {
+			delete dls ( i );
+		}
+	}
 };
 
@@ -430,6 +493,6 @@
 public:
 	//!Default constructor
-	mmix():iepdf(new emix()) {
-	    set_ep(iepdf);
+	mmix() : iepdf ( new emix() ) {
+		set_ep ( iepdf );
 	}
 
@@ -438,13 +501,15 @@
 		Array<epdf*> Eps ( Coms.length() );
 
-		for ( int i = 0;i < Coms.length();i++ ) {
-			Eps(i) = Coms(i)->e();
-		}
-
-		iepdf->set_parameters(w, Eps);
+		for ( int i = 0; i < Coms.length(); i++ ) {
+			Eps ( i ) = Coms ( i )->e();
+		}
+
+		iepdf->set_parameters ( w, Eps );
 	}
 
 	void condition ( const vec &cond ) {
-		for ( int i = 0;i < Coms.length();i++ ) {Coms ( i )->condition ( cond );}
+		for ( int i = 0; i < Coms.length(); i++ ) {
+			Coms ( i )->condition ( cond );
+		}
 	};
 };
Index: library/bdm/stat/exp_family.cpp
===================================================================
--- library/bdm/stat/exp_family.cpp (revision 471)
+++ library/bdm/stat/exp_family.cpp (revision 477)
@@ -4,5 +4,5 @@
 #include "exp_family.h"
 
-namespace bdm{
+namespace bdm {
 
 Uniform_RNG UniRNG;
@@ -12,5 +12,7 @@
 using std::cout;
 
-void BMEF::bayes ( const vec &dt ) {this->bayes ( dt,1.0 );};
+void BMEF::bayes ( const vec &dt ) {
+	this->bayes ( dt, 1.0 );
+};
 
 vec egiw::sample() const {
@@ -20,22 +22,21 @@
 
 double egiw::evallog_nn ( const vec &val ) const {
-	int vend = val.length()-1;
-
-	if ( dimx==1 ) { //same as the following, just quicker.
+	int vend = val.length() - 1;
+
+	if ( dimx == 1 ) { //same as the following, just quicker.
 		double r = val ( vend ); //last entry!
-		if (r<0) return -inf;
-		vec Psi ( nPsi+dimx );
+		if ( r < 0 ) return -inf;
+		vec Psi ( nPsi + dimx );
 		Psi ( 0 ) = -1.0;
-		Psi.set_subvector ( 1,val ( 0,vend-1 ) ); // fill the rest
-
-		double Vq=V.qform ( Psi );
-		return -0.5* ( nu*log ( r ) + Vq /r );
-	}
-	else {
-		mat Th= reshape ( val ( 0,nPsi*dimx-1 ),nPsi,dimx );
-		fsqmat R ( reshape ( val ( nPsi*dimx,vend ),dimx,dimx ) );
-		double ldetR=R.logdet();
-		if (ldetR) return -inf;
-		mat Tmp=concat_vertical ( -eye ( dimx ),Th );
+		Psi.set_subvector ( 1, val ( 0, vend - 1 ) ); // fill the rest
+
+		double Vq = V.qform ( Psi );
+		return -0.5* ( nu*log ( r ) + Vq / r );
+	} else {
+		mat Th = reshape ( val ( 0, nPsi * dimx - 1 ), nPsi, dimx );
+		fsqmat R ( reshape ( val ( nPsi*dimx, vend ), dimx, dimx ) );
+		double ldetR = R.logdet();
+		if ( ldetR ) return -inf;
+		mat Tmp = concat_vertical ( -eye ( dimx ), Th );
 		fsqmat iR ( dimx );
 		R.inv ( iR );
@@ -48,5 +49,5 @@
 	const vec& D = V._D();
 
-	double m = nu - nPsi -dimx-1;
+	double m = nu - nPsi - dimx - 1;
 #define	log2  0.693147180559945286226763983
 #define	logpi 1.144729885849400163877476189
@@ -54,29 +55,30 @@
 #define Inf std::numeric_limits<double>::infinity()
 
-	double nkG = 0.5* dimx* ( -nPsi *log2pi + sum ( log ( D ( dimx,D.length()-1 ) ) ) );
+	double nkG = 0.5 * dimx * ( -nPsi * log2pi + sum ( log ( D ( dimx, D.length() - 1 ) ) ) );
 	// temporary for lgamma in Wishart
-	double lg=0;
-	for ( int i =0; i<dimx;i++ ) {lg+=lgamma ( 0.5* ( m-i ) );}
-
-	double nkW = 0.5* ( m*sum ( log ( D ( 0,dimx-1 ) ) ) ) \
-	             - 0.5*dimx* ( m*log2 + 0.5* ( dimx-1 ) *log2pi )  - lg;
-
-	it_assert_debug ( ( ( -nkG-nkW ) >-Inf ) && ( ( -nkG-nkW ) <Inf ), "ARX improper" );
-	return -nkG-nkW;
+	double lg = 0;
+	for ( int i = 0; i < dimx; i++ ) {
+		lg += lgamma ( 0.5 * ( m - i ) );
+	}
+
+	double nkW = 0.5 * ( m * sum ( log ( D ( 0, dimx - 1 ) ) ) ) \
+	             - 0.5 * dimx * ( m * log2 + 0.5 * ( dimx - 1 ) * log2pi )  - lg;
+
+	it_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" );
+	return -nkG - nkW;
 }
 
 vec egiw::est_theta() const {
-	if ( dimx==1 ) {
+	if ( dimx == 1 ) {
 		const mat &L = V._L();
 		int end = L.rows() - 1;
 
-		mat iLsub = ltuinv ( L ( dimx,end,dimx,end ) );
-
-		vec L0 = L.get_col(0);
-
-		return iLsub * L0(1,end);
-	}
-	else {
-		it_error("ERROR: est_theta() not implemented for dimx>1");
+		mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) );
+
+		vec L0 = L.get_col ( 0 );
+
+		return iLsub * L0 ( 1, end );
+	} else {
+		it_error ( "ERROR: est_theta() not implemented for dimx>1" );
 		return 0;
 	}
@@ -89,12 +91,11 @@
 		int end = D.length() - 1;
 
-		mat Lsub = L ( 1, end, 1, end);
-		mat Dsub = diag(D(1, end));
-
-		return inv( transpose(Lsub) * Dsub * Lsub );
-
-	}
-	else {
-		it_error("ERROR: est_theta_cov() not implemented for dimx>1");
+		mat Lsub = L ( 1, end, 1, end );
+		mat Dsub = diag ( D ( 1, end ) );
+
+		return inv ( transpose ( Lsub ) * Dsub * Lsub );
+
+	} else {
+		it_error ( "ERROR: est_theta_cov() not implemented for dimx>1" );
 		return 0;
 	}
@@ -104,18 +105,17 @@
 vec egiw::mean() const {
 
-	if ( dimx==1 ) {
-		const vec &D= V._D();
-		int end = D.length()-1;
+	if ( dimx == 1 ) {
+		const vec &D = V._D();
+		int end = D.length() - 1;
 
 		vec m ( dim );
 		m.set_subvector ( 0, est_theta() );
-		m ( end ) = D ( 0 ) / ( nu -nPsi -2*dimx -2 );
+		m ( end ) = D ( 0 ) / ( nu - nPsi - 2 * dimx - 2 );
 		return m;
-	}
-	else {
+	} else {
 		mat M;
 		mat R;
-		mean_mat ( M,R );
-		return cvectorize ( concat_vertical ( M,R ) );
+		mean_mat ( M, R );
+		return cvectorize ( concat_vertical ( M, R ) );
 	}
 
@@ -124,43 +124,44 @@
 vec egiw::variance() const {
 
-	if ( dimx==1 ) {
-		int l=V.rows();
-		const ldmat tmp(V,linspace(1,l-1));
-		ldmat itmp(l);
-		tmp.inv(itmp);
-		double cove = V._D() ( 0 ) / ( nu -nPsi -2*dimx -2 );
-		
-		vec var(l);
-		var.set_subvector(0,diag(itmp.to_mat())*cove);
-		var(l-1)=cove*cove/( nu -nPsi -2*dimx -2 );
+	if ( dimx == 1 ) {
+		int l = V.rows();
+		const ldmat tmp ( V, linspace ( 1, l - 1 ) );
+		ldmat itmp ( l );
+		tmp.inv ( itmp );
+		double cove = V._D() ( 0 ) / ( nu - nPsi - 2 * dimx - 2 );
+
+		vec var ( l );
+		var.set_subvector ( 0, diag ( itmp.to_mat() ) *cove );
+		var ( l - 1 ) = cove * cove / ( nu - nPsi - 2 * dimx - 2 );
 		return var;
-	}
-	else {it_error("not implemented"); return vec(0);}
+	} else {
+		it_error ( "not implemented" );
+		return vec ( 0 );
+	}
 }
 
 void egiw::mean_mat ( mat &M, mat&R ) const {
-	const mat &L= V._L();
-	const vec &D= V._D();
-	int end = L.rows()-1;
-
-	ldmat ldR ( L ( 0,dimx-1,0,dimx-1 ), D ( 0,dimx-1 ) / ( nu -nPsi -2*dimx -2 ) ); //exp val of R
-	mat iLsub=ltuinv ( L ( dimx,end,dimx,end ) );
+	const mat &L = V._L();
+	const vec &D = V._D();
+	int end = L.rows() - 1;
+
+	ldmat ldR ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) / ( nu - nPsi - 2*dimx - 2 ) ); //exp val of R
+	mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) );
 
 	// set mean value
-	mat Lpsi = L ( dimx,end,0,dimx-1 );
-	M= iLsub*Lpsi;
-	R= ldR.to_mat()  ;
+	mat Lpsi = L ( dimx, end, 0, dimx - 1 );
+	M = iLsub * Lpsi;
+	R = ldR.to_mat()  ;
 }
 
 vec egamma::sample() const {
-	vec smp ( dim);
+	vec smp ( dim );
 	int i;
 
-	for ( i=0; i<dim; i++ ) {
-		if ( beta ( i ) >std::numeric_limits<double>::epsilon() ) {
-			GamRNG.setup ( alpha ( i ),beta ( i ) );
-		}
-		else {
-			GamRNG.setup ( alpha ( i ),std::numeric_limits<double>::epsilon() );
+	for ( i = 0; i < dim; i++ ) {
+		if ( beta ( i ) > std::numeric_limits<double>::epsilon() ) {
+			GamRNG.setup ( alpha ( i ), beta ( i ) );
+		} else {
+			GamRNG.setup ( alpha ( i ), std::numeric_limits<double>::epsilon() );
 		}
 #pragma omp critical
@@ -190,11 +191,11 @@
 	int i;
 
-	if (any(val<=0.)) return -inf;
-	if (any(beta<=0.)) return -inf;
-	for ( i=0; i<dim; i++ ) {
-		res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) - beta ( i ) *val ( i );
-	}
-	double tmp=res-lognc();;
-	it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+	if ( any ( val <= 0. ) ) return -inf;
+	if ( any ( beta <= 0. ) ) return -inf;
+	for ( i = 0; i < dim; i++ ) {
+		res += ( alpha ( i ) - 1 ) * std::log ( val ( i ) ) - beta ( i ) * val ( i );
+	}
+	double tmp = res - lognc();;
+	it_assert_debug ( std::isfinite ( tmp ), "Infinite value" );
 	return tmp;
 }
@@ -204,6 +205,6 @@
 	int i;
 
-	for ( i=0; i<dim; i++ ) {
-		res += lgamma ( alpha ( i ) ) - alpha ( i ) *std::log ( beta ( i ) ) ;
+	for ( i = 0; i < dim; i++ ) {
+		res += lgamma ( alpha ( i ) ) - alpha ( i ) * std::log ( beta ( i ) ) ;
 	}
 
@@ -211,47 +212,47 @@
 }
 
-void mgamma::set_parameters(double k0, const vec &beta0) {
+void mgamma::set_parameters ( double k0, const vec &beta0 ) {
 	k = k0;
-	iepdf->set_parameters(k * ones(beta0.length()), beta0);
+	iepdf->set_parameters ( k * ones ( beta0.length() ), beta0 );
 	dimc = e()->dimension();
 }
 
-ivec eEmp::resample (RESAMPLING_METHOD method) {
-	ivec ind=zeros_i ( n );
+ivec eEmp::resample ( RESAMPLING_METHOD method ) {
+	ivec ind = zeros_i ( n );
 	ivec N_babies = zeros_i ( n );
 	vec cumDist = cumsum ( w );
 	vec u ( n );
-	int i,j,parent;
+	int i, j, parent;
 	double u0;
 
 	switch ( method ) {
-		case MULTINOMIAL:
-			u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
-
-			for ( i = n - 2;i >= 0;i-- ) {
-				u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
-			}
-
-			break;
-
-		case STRATIFIED:
-
-			for ( i = 0;i < n;i++ ) {
-				u ( i ) = ( i + UniRNG.sample() ) / n;
-			}
-
-			break;
-
-		case SYSTEMATIC:
-			u0 = UniRNG.sample();
-
-			for ( i = 0;i < n;i++ ) {
-				u ( i ) = ( i + u0 ) / n;
-			}
-
-			break;
-
-		default:
-			it_error ( "PF::resample(): Unknown resampling method" );
+	case MULTINOMIAL:
+		u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
+
+		for ( i = n - 2; i >= 0; i-- ) {
+			u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
+		}
+
+		break;
+
+	case STRATIFIED:
+
+		for ( i = 0; i < n; i++ ) {
+			u ( i ) = ( i + UniRNG.sample() ) / n;
+		}
+
+		break;
+
+	case SYSTEMATIC:
+		u0 = UniRNG.sample();
+
+		for ( i = 0; i < n; i++ ) {
+			u ( i ) = ( i + u0 ) / n;
+		}
+
+		break;
+
+	default:
+		it_error ( "PF::resample(): Unknown resampling method" );
 	}
 
@@ -259,5 +260,5 @@
 	j = 0;
 
-	for ( i = 0;i < n;i++ ) {
+	for ( i = 0; i < n; i++ ) {
 		while ( u ( i ) > cumDist ( j ) ) j++;
 
@@ -270,16 +271,16 @@
 
 	// find the first parent;
-	parent=0; while ( N_babies ( parent ) ==0 ) parent++;
+	parent = 0;
+	while ( N_babies ( parent ) == 0 ) parent++;
 
 	// Build index
-	for ( i = 0;i < n;i++ ) {
+	for ( i = 0; i < n; i++ ) {
 		if ( N_babies ( i ) > 0 ) {
 			ind ( i ) = i;
 			N_babies ( i ) --; //this index was now replicated;
-		}
-		else {
+		} else {
 			// test if the parent has been fully replicated
 			// if yes, find the next one
-			while ( ( N_babies ( parent ) ==0 ) || ( N_babies ( parent ) ==1 && parent>i ) ) parent++;
+			while ( ( N_babies ( parent ) == 0 ) || ( N_babies ( parent ) == 1 && parent > i ) ) parent++;
 
 			// Replicate parent
@@ -292,9 +293,9 @@
 
 	// copy the internals according to ind
-	for ( i=0;i<n;i++ ) {
-		if ( ind ( i ) !=i ) {
-			samples ( i ) =samples ( ind ( i ) );
-		}
-		w ( i ) = 1.0/n;
+	for ( i = 0; i < n; i++ ) {
+		if ( ind ( i ) != i ) {
+			samples ( i ) = samples ( ind ( i ) );
+		}
+		w ( i ) = 1.0 / n;
 	}
 
@@ -305,34 +306,36 @@
 	//it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
 	dim = epdf0->dimension();
-	w=w0;
-	w/=sum ( w0 );//renormalize
-	n=w.length();
+	w = w0;
+	w /= sum ( w0 );//renormalize
+	n = w.length();
 	samples.set_size ( n );
 	dim = epdf0->dimension();
 
-	for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}
+	for ( int i = 0; i < n; i++ ) {
+		samples ( i ) = epdf0->sample();
+	}
 }
 
 void eEmp::set_samples ( const epdf* epdf0 ) {
 	//it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
-	w=1;
-	w/=sum ( w );//renormalize
-
-	for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}
-}
-
-void migamma_ref::from_setting( const Setting &set ) 
-{				
+	w = 1;
+	w /= sum ( w );//renormalize
+
+	for ( int i = 0; i < n; i++ ) {
+		samples ( i ) = epdf0->sample();
+	}
+}
+
+void migamma_ref::from_setting ( const Setting &set ) {
 	vec ref;
-	UI::get( ref, set, "ref" , UI::compulsory);
-	set_parameters(set["k"],ref,set["l"]);
-}
-
-void mlognorm::from_setting( const Setting &set ) 
-{	
-	vec mu0;	
-	UI::get( mu0, set, "mu0", UI::compulsory);
-	set_parameters(mu0.length(),set["k"]);
-	condition(mu0);
+	UI::get ( ref, set, "ref" , UI::compulsory );
+	set_parameters ( set["k"], ref, set["l"] );
+}
+
+void mlognorm::from_setting ( const Setting &set ) {
+	vec mu0;
+	UI::get ( mu0, set, "mu0", UI::compulsory );
+	set_parameters ( mu0.length(), set["k"] );
+	condition ( mu0 );
 }
 
