Index: /bdm/estim/merger.cpp
===================================================================
--- /bdm/estim/merger.cpp (revision 311)
+++ /bdm/estim/merger.cpp (revision 379)
@@ -5,25 +5,36 @@
 namespace bdm
 {
-	vec merger::lognorm_merge ( mat &lW )
+	vec merger_base::merge_points ( mat &lW )
 	{
 		int nu=lW.rows();
-		vec mu = sum ( lW ) /nu; //mean of logs
-		// vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); ======= numerically unsafe!
-		vec 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 )
+		switch ( METHOD )
 		{
-			case 2:
-				coef= ( 1-0.5*sqrt ( ( 4.0*beta-3.0 ) /beta ) );
-				return  coef*sq2bl + mu ;
+			case ARITHMETIC:
+				return log ( sum ( exp ( lW ) ) ); //ugly!
 				break;
-			case 3://Ratio of Bessel
-				coef = sqrt ( ( 3*beta-2 ) /3*beta );
-				return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) +  mu;
+			case GEOMETRIC:
+				return sum ( lW ) /nu;
 				break;
-			case 4:
-				break;
-			default: // Approximate conditional density
+			case LOGNORMAL:
+				vec mu = sum ( lW ) /nu; //mean of logs
+				// vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); ======= numerically unsafe!
+				vec 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 ) );
+						return  coef*sq2bl + mu ;
+						break;
+					case 3://Ratio of Bessel
+						coef = sqrt ( ( 3*beta-2 ) /3*beta );
+						return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) +  mu;
+						break;
+					case 4:
+						break;
+					default: // Approximate conditional density
+						break;
+				}
 				break;
 		}
@@ -31,35 +42,23 @@
 	}
 
-	void merger::merge ( const epdf* g0 )
+	void merger_mix::merge ( )
 	{
-
-		it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" );
-		//Empirical density - samples
-		if ( !fix_smp )
-		{
-			eSmp.set_statistics ( ones ( Ns ), g0 );
-		}
-
 		Array<vec> &Smp = eSmp._samples(); //aux
 		vec &w = eSmp._w(); //aux
 
-		mat Smp_ex =ones ( dim +1,Ns ); // Extended samples for the ARX model - the last row is ones
-		for ( int i=0;i<Ns;i++ ) {	set_col_part ( Smp_ex,i,Smp ( i ) );}
-
-		if ( DBG )	*dbg << Name ( "Smp_0" ) << Smp_ex;
+		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 ( Ns );
-		vec lw_mix ( Ns );
-		vec lw ( Ns );
-		mat lW=zeros ( n,Ns );
+		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
-		if ( !fix_smp )
-			for ( int i=0;i<Ns;i++ )
-			{
-				lw_mix ( i ) =g0->evallog ( Smp ( i ) );
-			}
+		lw_mix = 1.0; // assuming uniform grid density -- otherwise 
 
 		// Initial component in the mixture model
@@ -67,5 +66,5 @@
 		ARX A0;	A0.set_statistics ( dim, V0 ); //initial guess of Mix:
 
-		Mix.init ( &A0, Smp_ex, Nc );
+		Mix.init ( &A0, Smp_ex, Ncoms );
 		//Preserve initial mixture for repetitive estimation via flattening
 		MixEF Mix_init ( Mix );
@@ -74,8 +73,8 @@
 		bool converged=false;
 		int niter = 0;
-		char str[100];
+		char dbg_str[100];
 
 		emix* Mpred=Mix.epredictor ( );
-		vec Mix_pdf ( Ns );
+		vec Mix_pdf ( Npoints );
 		while ( !converged )
 		{
@@ -83,5 +82,5 @@
 			//Re-Initialize Mixture model
 			Mix.flatten ( &Mix_init );
-			Mix.bayesB ( Smp_ex, w*Ns );
+			Mix.bayesB ( Smp_ex, w*Npoints );
 			delete Mpred;
 			Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!!
@@ -89,12 +88,12 @@
 
 			// This will be active only later in iterations!!!
-			if ( ( !fix_smp ) & ( 1./sum_sqr ( w ) <effss_coef*Ns ) )
+			if (  1./sum_sqr ( w ) <effss_coef*Npoints ) 
 			{
 				// Generate new samples
 				eSmp.set_samples ( Mpred );
-				for ( int i=0;i<Ns;i++ )
+				for ( int i=0;i<Npoints;i++ )
 				{
 					//////////// !!!!!!!!!!!!!
-					if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; }
+					//if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; }
 					set_col_part ( Smp_ex,i,Smp ( i ) );
 					//Importance of the mixture
@@ -102,30 +101,30 @@
 					lw_mix ( i ) = Mpred->evallog ( Smp ( i ) );
 				}
-				if ( 0 )
+				if ( DBG )
 				{
 					cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;
 					cout << Mix._e()->mean() <<endl;
-					cout << sum ( Smp_ex,2 ) /Ns <<endl;
-					cout << Smp_ex*Smp_ex.T() /Ns << endl;
+					cout << sum ( Smp_ex,2 ) /Npoints <<endl;
+					cout << Smp_ex*Smp_ex.T() /Npoints << endl;
 				}
 			}
 			if ( DBG )
 			{
-				sprintf ( str,"Mpred_mean%d",niter );
-				*dbg << Name ( str ) << Mpred->mean();
-				sprintf ( str,"Mpred_var%d",niter );
-				*dbg << Name ( str ) << Mpred->variance();
-
-
-				sprintf ( str,"Mpdf%d",niter );
-				for ( int i=0;i<Ns;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );}
-				*dbg << Name ( str ) << Mix_pdf;
-
-				sprintf ( str,"Smp%d",niter );
-				*dbg << Name ( str ) << Smp_ex;
+				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<n;i++ )
+			for ( int i=0;i<mpdfs.length();i++ )
 			{
 				lw_src=0.0;
@@ -135,5 +134,5 @@
 				{
 					// no need for conditioning or marginalization
-					for ( int j=0;j<Ns; j++ )   // Smp is Array<> => for cycle
+					for ( int j=0;j<Npoints; j++ )   // Smp is Array<> => for cycle
 					{
 						lw_src ( j ) =mpdfs ( i )->_epdf().evallog ( Smp ( j ) );
@@ -148,5 +147,5 @@
 						epdf* tmp_marg = Mpred->marginal ( mpdfs ( i )->_rvc() );
 						//compute vector of lw_src
-						for ( int k=0;k<Ns;k++ )
+						for ( int k=0;k<Npoints;k++ )
 						{
 							// Here val of tmp_marg = cond of mpdfs(i) ==> calling dls->get_cond
@@ -167,5 +166,5 @@
 						// Compute likelihood
 						vec lw_dbg=lw_src;
-						for ( int k= 0; k<Ns; k++ )
+						for ( int k= 0; k<Npoints; k++ )
 						{
 							lw_src ( k ) += log (
@@ -181,5 +180,5 @@
 					}
 					// Compute likelihood of the partial source
-					for ( int k= 0; k<Ns; k++ )
+					for ( int k= 0; k<Npoints; k++ )
 					{
 						mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp ( k ) ) );
@@ -188,12 +187,11 @@
 
 				}
-//			it_assert_debug(std::isfinite(sum(lw_src)),"bad");
+		//			it_assert_debug(std::isfinite(sum(lw_src)),"bad");
 				lW.set_row ( i, lw_src ); // do not divide by mix
 			}
-			lw = lognorm_merge ( lW ); //merge
+			lw = merger_base::merge_points ( lW ); //merge
 
 			//Importance weighting
-			if ( !fix_smp )
-				lw -=  lw_mix; // hoping that it is not numerically sensitive...
+			lw -=  lw_mix; // hoping that it is not numerically sensitive...
 			w = exp ( lw-max ( lw ) );
 
@@ -212,10 +210,10 @@
 			if ( DBG )
 			{
-				sprintf ( str,"lW%d",niter );
-				*dbg << Name ( str ) << lW;
-				sprintf ( str,"w%d",niter );
-				*dbg << Name ( str ) << w;
-				sprintf ( str,"lw_m%d",niter );
-				*dbg << Name ( str ) << lw_mix;
+				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 ===
Index: /bdm/estim/merger.h
===================================================================
--- /bdm/estim/merger.h (revision 311)
+++ /bdm/estim/merger.h (revision 379)
@@ -19,190 +19,257 @@
 namespace bdm
 {
-	using std::string;
-
-	/*!
-	@brief Function for general combination of pdfs
-
-	Mixtures of Gaussian densities are used internally. Switching to other densities should be trivial.
-	*/
-
-	class merger : public compositepdf, public epdf
-	{
-		protected:
-			//!Internal mixture of EF models
-			MixEF Mix;
-			//! Data link for each mpdf in mpdfs
-			Array<datalink_m2e*> dls;
-			//! Array of rvs that are not modelled by mpdfs at all (aux)
-			Array<RV> rvzs;
-			//! Data Links of rv0 mpdfs - these will be conditioned the (rv,rvc) of mpdfs
-			Array<datalink_m2e*> zdls;
-
-			//!Number of samples used in approximation
-			int Ns;
-			//!Number of components in a mixture
-			int Nc;
-			//!Prior on the log-normal merging model
-			double beta;
-			//! Projection to empirical density
-			eEmp eSmp;
-			//! coefficient of resampling
-			double effss_coef;
-
-			//! debug or not debug
-			bool DBG;
-			//! debugging file
-			it_file* dbg;
-			//! Flag if the samples are fixed or not
-			bool fix_smp;
-		public:
-//!Default constructor
-			merger ( const Array<mpdf*> &S ) :
-					compositepdf ( S ), epdf ( ),
-					Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ), eSmp()
-			{
-				RV ztmp;
-				rv = getrv ( false );
-				RV rvc; setrvc ( rv,rvc ); // Extend rv by rvc!
-				rv.add ( rvc );
-				// get dimension
-				dim = rv._dsize();
-
-				for ( int i=0;i<n;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
-					ztmp= mpdfs ( i )->_rv();
-					ztmp.add ( mpdfs ( i )->_rvc() );
-					rvzs ( i ) =rv.subt ( ztmp );
-					zdls ( i ) = new datalink_m2e; zdls ( i )->set_connection ( rvzs ( i ), ztmp, rv ) ;
-				};
-				//Set Default values of parameters
-				beta=2.0;
-				Ns=100;
-				Nc=10;
-				Mix.set_method ( EM );
-				DBG = false;
-				fix_smp = false;
-			}
-			//! set debug file
-			void debug_file ( const string fname ) { if ( DBG ) delete dbg; dbg = new it_file ( fname ); if ( dbg ) DBG=true;}
-//! Set internal parameters used in approximation
-			void set_parameters ( double beta0, int Ns0, int Nc0, double effss_coef0=0.5 ) {beta=beta0;
-				Ns=Ns0;
-				Nc=Nc0;
-				effss_coef=effss_coef0;
-				eSmp.set_parameters ( Ns0,false );
-			}
-			void set_grid ( Array<vec> &XYZ )
-			{
-				int dim=XYZ.length(); ivec szs ( dim );
-				for(int i=0; i<dim;i++){szs=XYZ(i).length();}
-				Ns=prod(szs);
-				eSmp.set_parameters(Ns,false);
-				Array<vec> &samples=eSmp._samples();
-				eSmp._w()=ones(Ns)/Ns;
-						
-				//set samples
-				ivec is=zeros_i(dim);//indeces of dimensions in for cycle;
-				vec smpi(dim);
-				for(int i=0; i<Ns; i++){
-					for(int j=0; j<dim; j++){smpi(j)=XYZ(j)(is(j)); /* jty vector*/ }
-					samples(i)=smpi;
-					// shift indeces
-					for (int j=0;j<dim;j++){
-						if (is(j)==szs(j)-1) { //j-th index is full
-							is(j)=0; //shift back
-							is(j+1)++; //increase th next dimension;
-							if (is(j+1)<szs(j+1)-1) break;
-						} else {
-							is(j)++; break;
-						}
+using std::string;
+
+//!Merging methods
+enum MERGER_METHOD {ARITHMETIC = 1, GEOMETRIC = 2, LOGNORMAL = 3};
+
+/*!
+@brief Base class for general combination of pdfs on discrete support
+
+Mixtures of Gaussian densities are used internally. Switching to other densities should be trivial.
+
+The merged pdfs are expected to be of the form:
+ \f[ f(x_i|y_i),  i=1..n \f]
+where the resulting merger is a density on \f$ \cup [x_i,y_i] \f$ .
+Note that all variables will be joined.
+
+As a result of this feature, each source must be extended to common support
+\f[ f(z_i|y_i,x_i) f(x_i|y_i) f(y_i)  i=1..n \f]
+where \f$ z_i \f$ accumulate variables that were not in the original source.
+These extensions are calculated on-the-fly.
+
+However, these operations can not be performed in general. Hence, this class merges only sources on common support, \f$ y_i={}, z_i={}, \forall i \f$.
+For merging of more general cases, use offsprings merger_mix and merger_grid.
+*/
+
+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;
+		//!Prior on the log-normal merging model
+		double 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) {set_sources (S);};
+		//! Function setting the main internal structures 
+		void set_sources (const Array<mpdf*> &Sources) {
+			compositepdf::set_elements (Sources);
+			//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) ;
+			};
+		}
+		//! 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, double beta0=0.0) {
+			METHOD = MTH; 
+			beta = beta0;
+		}
+		//! Set support points from a pdf by drawing N samples
+		void set_support(const epdf &overall, int N){
+			it_assert_debug ( rv.equal ( overall._rv() ),"Incompatible parameter overall!" );
+			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
+		void merge () {
+			if (eSmp._w().length() ==0) {it_error("Empty support points use set_support" );}
+			//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);
 					}
 				}
-				
-				fix_smp = true;
-			}
-//!Initialize the proposal density. This function must be called before merge()!
-			void init()   ////////////// NOT FINISHED
-			{
-				Array<vec> Smps ( n );
-				//Gibbs sampling
-				for ( int i=0;i<n;i++ ) {Smps ( i ) =zeros ( 0 );}
-			}
-//!Create a mixture density using known proposal
-			void merge ( const epdf* g0 );
-//!Create a mixture density, make sure to call init() before the first call
-			void merge () {merge ( & ( Mix.posterior() ) );};
-
-//! Merge log-likelihood values
-			vec lognorm_merge ( mat &lW );
-//! sample from merged density
+
+				vec wtmp = exp (merge_points (lW));
+				//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 sample ( ) const { return Mix.posterior().sample();}
-			double evallog ( const vec &dt ) const
-			{
-				vec dtf=ones ( dt.length() +1 );
-				dtf.set_subvector ( 0,dt );
-				return Mix.logpred ( dtf );
-			}
-			vec mean() const
-			{
-				const Vec<double> &w = eSmp._w();
-				const Array<vec> &S = eSmp._samples();
-				vec tmp=zeros ( dim );
-				for ( int i=0; i<Ns; 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<Ns; 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<Ns; i++ )
-				{
-					tmp+=w ( i ) *pow ( S ( i ),2 );
-				}
-				return tmp-pow ( mean(),2 );
-			}
-//! for future use
-			virtual ~merger()
-			{
-				for ( int i=0; i<n; i++ )
-				{
-					delete dls ( i );
-					delete zdls ( i );
-				}
-				if ( DBG ) delete dbg;
-			};
-
+		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;}
+		
+		//!@}
+};
+
+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
+		double effss_coef;
+
+	public:
+		//!\name Constructors
+		//!@{
+		merger_mix () {};
+		merger_mix (const Array<mpdf*> &S) {set_sources(S);};
+		//! Set sources and prepare all internal structures
+		void set_sources (const Array<mpdf*> &S) {
+			merger_base::set_sources(S);
+			Nsources = S.length();
+		}
+		//! Set internal parameters used in approximation
+		void set_parameters (int Ncoms0=10, double effss_coef0 = 0.5) {
+			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;}
-//! Access function
-			eEmp& _Smp() {return eSmp;}
-	};
+		MixEF& _Mix() {return Mix;}
+		//! Access function
+		emix* proposal() {emix* tmp = Mix.epredictor(); tmp->set_rv (rv); return tmp;}
+	//! @}
+		
+};
 
 }
Index: /applications/mpdm/TR2244/merger_iter_cond.cpp
===================================================================
--- /applications/mpdm/TR2244/merger_iter_cond.cpp (revision 310)
+++ /applications/mpdm/TR2244/merger_iter_cond.cpp (revision 379)
@@ -34,6 +34,6 @@
 	A ( 1 ) =&mf2;
 
-	merger M ( A );
-	M.debug_file("iter_cond_debug.it");
+	merger_mix M ( A );
+	M.set_debug_file("iter_cond_debug.it");
 	enorm<ldmat> g0; g0.set_rv ( all );
 	mat Cov=10*eye ( 2 );
@@ -45,5 +45,6 @@
 //	mlnorm<ldmat>* testm=(mlnorm<ldmat>*)teste->condition(u2);
 	
-	M.set_parameters ( 1.2,1000,1,0.5);
+	M.set_parameters (1,0.5);
+	M.set_method(LOGNORMAL, 1.2);
 	
 	Array<vec> YUU(2); 
@@ -66,6 +67,7 @@
 
 	for ( int it=0;it<Ntrials;it++ ) {
-		M.merge ( &g0 );
-
+		M.set_support( g0,1000);
+		M.merge();
+		
 		MixEF &MM = M._Mix();
 		epdf* MP = MM._Coms ( 0 )->epredictor ( );
Index: /applications/mpdm/SYSID09/merg_2a.cpp
===================================================================
--- /applications/mpdm/SYSID09/merg_2a.cpp (revision 321)
+++ /applications/mpdm/SYSID09/merg_2a.cpp (revision 379)
@@ -126,12 +126,14 @@
 		mepdf eG2 ( P2._e() );
 		Array<mpdf*> A ( 2 ); A ( 0 ) =&eG1;A ( 1 ) =&eG2;
-		merger M ( A );
-		M.set_parameters ( 2, 100,20 ,0.99); M._Mix().set_method(QB);
+		merger_mix M ( A );
+		M.set_parameters ( 20 ,0.99); 
+		M.set_method(LOGNORMAL, 1.2);
+		M._Mix().set_method(QB);
 		//M2.set_parameters ( 100.0, 1000,3 ); //M2._Mix().set_method(QB);
 /*		char fnm[100];
 		sprintf(fnm,"m2a_dbg%d.it",t);
 		M.debug_file(fnm);*/
-		M.merge ( proposal );
-
+		M.set_support ( *proposal,100 );
+		M.merge();
 		//proposal = M.proposal();
 		//Likelihood
