Index: /bdm/estim/merger.cpp
===================================================================
--- /bdm/estim/merger.cpp (revision 211)
+++ /bdm/estim/merger.cpp (revision 213)
@@ -6,5 +6,6 @@
 	int nu=lW.rows();
 	vec mu = sum ( lW ) /nu; //mean of logs
-	vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 );
+	// 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
@@ -61,5 +62,5 @@
 	char str[100];
 
-	epdf* Mpred=Mix.predictor(rv);
+	epdf* Mpred=Mix.predictor ( rv );
 	vec Mix_pdf ( Ns );
 	while ( !converged ) {
@@ -70,5 +71,6 @@
 		delete Mpred;
 		Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!!
-		
+
+		// This will be active only later in iterations!!!
 		if ( 1./sum_sqr ( w ) <0.5*Ns ) {
 			// Generate new samples
@@ -76,10 +78,10 @@
 			for ( int i=0;i<Ns;i++ ) {
 				//////////// !!!!!!!!!!!!!
-				if ( Smp ( i ) ( 1 ) <0 ) {Smp ( i ) ( 1 ) = 0.01; }
+				if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; }
 				set_col_part ( Smp_ex,i,Smp ( i ) );
 			}
-			{cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;}
+			if(0){cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;
 			cout << sum ( Smp_ex,2 ) /Ns <<endl;
-			cout << Smp_ex*Smp_ex.T()/Ns << endl; 
+			cout << Smp_ex*Smp_ex.T() /Ns << endl;}
 		}
 // 		sprintf ( str,"Mpred_mean%d",niter );
@@ -172,5 +174,6 @@
 		}
 		else {
-
+			it_file itf ( "merg_err.it" );
+			itf << Name ( "w" ) << w;
 		}
 
@@ -188,4 +191,5 @@
 		converged = ( niter>20 );
 	}
+	delete Mpred;
 	cout << endl;
 
Index: /tests/merger_2d_test.cpp
===================================================================
--- /tests/merger_2d_test.cpp (revision 211)
+++ /tests/merger_2d_test.cpp (revision 213)
@@ -44,5 +44,5 @@
 	merger M ( A );
 	enorm<fsqmat> g0(xy);
-	g0.set_parameters(vec("1 1"),mat("1 0; 0 1"));
+	g0.set_parameters(vec("1 1"),mat("10 0; 0 10"));
 	
 	M.set_parameters(1e8,400,1);
Index: /CMakeLists.txt
===================================================================
--- /CMakeLists.txt (revision 211)
+++ /CMakeLists.txt (revision 213)
@@ -27,5 +27,6 @@
 	#This is for UNIX makefile which does only one release at a time.
 	SET(CMAKE_BUILD_TYPE Debug)
-	SET(CMAKE_CXX_FLAGS_DEBUG "-Wall -g -DNDEBUG -O3 -pipe  -Wno-unknown-pragmas")
+#	SET(CMAKE_CXX_FLAGS_DEBUG "-Wall -g -DNDEBUG -O3 -pipe  -Wno-unknown-pragmas")
+	SET(CMAKE_CXX_FLAGS_DEBUG "-Wall -g -pipe  -Wno-unknown-pragmas")
 	INCLUDE(CMakeLocal.txt OPTIONAL)	
 ENDIF(WIN32)
Index: /mpdm/merg_pred.cpp
===================================================================
--- /mpdm/merg_pred.cpp (revision 211)
+++ /mpdm/merg_pred.cpp (revision 213)
@@ -15,7 +15,7 @@
 	RV u1 ( "{u1 }" );
 	RV u2 ( "{u2 }" );
-	RV ym=y; ym.t(-1);
-	RV yy = y; yy.add(ym);
-	RV uu=u1; uu.add(u2);
+	RV ym=y; ym.t ( -1 );
+	RV yy = y; yy.add ( ym );
+	RV uu=u1; uu.add ( u2 );
 
 	// Full system
@@ -44,5 +44,5 @@
 
 	// Logging
-	dirfilelog L ( "exp/merg",ndat );
+	dirfilelog L ( "exp/merg",10 );
 
 	int Li_Eth1 = L.add ( thri,"P1" );
@@ -56,14 +56,14 @@
 
 	vec Yt ( ndat );
-	vec yt(1);
+	vec yt ( 1 );
 
 	Yt.set_subvector ( 0,randn ( ord ) ); //initial values
-	vec rgrg ( thrg.count() -1); // constant terms are in!
-	vec rgr1 ( thri.count() -1);
-	vec rgr2 ( thri.count() -1);
+	vec rgrg ( thrg.count() -1 ); // constant terms are in!
+	vec rgr1 ( thri.count() -1 );
+	vec rgr2 ( thri.count() -1 );
 
-	vec PredLLs(5);
-	vec PostLLs(3);
-	vec PredLLs_m=zeros(5);
+	vec PredLLs ( 5 );
+	vec PostLLs ( 3 );
+	vec PredLLs_m=zeros ( 5 );
 	ivec ind_r1 = "0 1 3";
 	ivec ind_r2 = "0 2 3";
@@ -73,31 +73,31 @@
 			rgrg ( 0 ) =Yt ( t-1 );
 			rgrg ( 1 ) = pow(sin ( ( t/40.0 ) *pi ),3);
-			rgrg ( 2 ) = pow(cos ( ( t/40.0 ) *pi ),3);
-			rgrg (3) = 1.0; // constant term
-			
-			rgr1(0) = rgrg(0); rgr1(1) = rgrg(1); rgr1(2) = rgrg(3); // no u2
-			rgr2(0) = rgrg(0); rgr2(1) = rgrg(2); rgr2(2) = rgrg(3); // no u1
+			rgrg ( 2 ) = pow(cos ( ( t/60.0 ) *pi ),3);
+			rgrg ( 3 ) = 1.0; // constant term
+
+			rgr1 ( 0 ) = rgrg ( 0 ); rgr1 ( 1 ) = rgrg ( 1 ); rgr1 ( 2 ) = rgrg ( 3 ); // no u2
+			rgr2 ( 0 ) = rgrg ( 0 ); rgr2 ( 1 ) = rgrg ( 2 ); rgr2 ( 2 ) = rgrg ( 3 ); // no u1
 
 			Yt ( t ) = thg*rgrg + sqr * NorRNG();
 
 			// Test predictors
-			if (t>2){
-				mlnorm<ldmat>* P1p = P1.predictor(y,concat(ym,u1));
-				mlnorm<ldmat>* P2p = P2.predictor(y,concat(ym,u2));
-				mlnorm<ldmat>* Pgp = PG.predictor(y,concat(ym,uu));
-				
-				Array<mpdf*> A(2); A(0)=P1p;A(1)=P2p;
-				merger M(A);
-				enorm<ldmat> g0(concat(yy,uu)); g0.set_parameters("0 0 0 0 ",3*eye(4));
-				M.set_parameters(10000000.0, 100,1);
-				M.merge(&g0);
-				
-				yt(0) = Yt(t);
-				double P1pl = P1p->evallogcond(yt,rgr1);	
-				double P2pl = P2p->evallogcond(yt,rgr2);	
-				double PGpl = Pgp->evallogcond(yt,rgrg);	
+			if ( t>2 ) {
+					mlnorm<ldmat>* P1p = P1.predictor ( y,concat ( ym,u1 ) );
+					mlnorm<ldmat>* P2p = P2.predictor ( y,concat ( ym,u2 ) );
+					mlnorm<ldmat>* Pgp = PG.predictor ( y,concat ( ym,uu ) );
+
+					Array<mpdf*> A ( 2 ); A ( 0 ) =P1p;A ( 1 ) =P2p;
+					merger M ( A );
+					enorm<ldmat> g0 ( concat ( yy,uu ) ); g0.set_parameters ( "0 0 0 0 ",3*eye ( 4 ) );
+					M.set_parameters ( 1e8, 101,1 );
+					M.merge ( &g0 );
+				yt ( 0 ) = Yt ( t );
+				double P1pl = P1p->evallogcond ( yt,rgr1 );
+				double P2pl = P2p->evallogcond ( yt,rgr2 );
+				double PGpl = Pgp->evallogcond ( yt,rgrg );
 				{
-					cout << "yt: " << yt << endl;
+					cout << "T: " << t << "yt: " << yt << endl;
 					cout << "yt_1: " << P1p->_epdf().mean() << endl;
+					cout << *P1p <<endl;
 					cout << "yt_2: " << P2p->_epdf().mean() << endl;
 					cout << "yt_G: " << P2p->_epdf().mean() << endl;
@@ -106,35 +106,37 @@
 				double cP2pl;
 				{
-					ARX* Apred = (ARX*)M._Mix()._Coms(0);
-					enorm<ldmat>* MP= Apred->predictor(concat(yy,uu));
-					enorm<ldmat>* mP1p = (enorm<ldmat>*)MP->marginal(concat(yy,u1));
-					enorm<ldmat>* mP2p = (enorm<ldmat>*)MP->marginal(concat(yy,u2));
-					mlnorm<ldmat>* cP1p = (mlnorm<ldmat>*)mP1p->condition(y);
-					mlnorm<ldmat>* cP2p = (mlnorm<ldmat>*)mP2p->condition(y);
+					ARX* Apred = ( ARX* ) M._Mix()._Coms ( 0 );
+					enorm<ldmat>* MP= Apred->predictor ( concat ( yy,uu ) );
+					enorm<ldmat>* mP1p = ( enorm<ldmat>* ) MP->marginal ( concat ( yy,u1 ) );
+					enorm<ldmat>* mP2p = ( enorm<ldmat>* ) MP->marginal ( concat ( yy,u2 ) );
+					mlnorm<ldmat>* cP1p = ( mlnorm<ldmat>* ) mP1p->condition ( y );
+					mlnorm<ldmat>* cP2p = ( mlnorm<ldmat>* ) mP2p->condition ( y );
 
-					cP1pl = cP1p->evallogcond(yt,rgr1);	
-					cP2pl = cP2p->evallogcond(yt,rgr2);	
-				
+					cP1pl = cP1p->evallogcond ( yt,rgr1 );
+					cP2pl = cP2p->evallogcond ( yt,rgr2 );
+
 					cout << "ytm1: " << cP1p->_epdf().mean() << endl;
 					cout << "ytm2: " << cP2p->_epdf().mean() << endl;
+					cout << *cP2p << endl;
 				}
 
 				PredLLs *=frg;
-				PredLLs += log(concat(vec_3(P1pl, P2pl, PGpl), vec_2(cP1pl,cP2pl)));
-				L.logit(Li_Pred, PredLLs); //log-normal
-				
+				PredLLs += concat ( vec_3 ( P1pl, P2pl, PGpl ), vec_2 ( cP1pl,cP2pl ) );
+				L.logit ( Li_Pred, PredLLs ); //log-normal
+
 				delete P1p;
 				delete P2p;
 				delete Pgp;
 			}
-			
-			// 1st
-			P1.bayes ( concat(Yt(t),rgr1) );
-			// 2nd
-			P2.bayes ( concat(Yt(t),rgr2) );
 
-			//Global
-			PG.bayes ( concat ( Yt ( t ),rgrg ) );
-			
+			if ( t<100 ) {
+				// 1st
+				P1.bayes ( concat ( Yt ( t ),rgr1 ) );
+				// 2nd
+				P2.bayes ( concat ( Yt ( t ),rgr2 ) );
+
+				//Global
+				PG.bayes ( concat ( Yt ( t ),rgrg ) );
+			}
 			//Merger
 		}
@@ -146,5 +148,5 @@
 		PostLLs += vec_3 ( P1._ll(), P2._ll(), PG._ll() );
 		L.logit ( Li_LL, PostLLs );
-		L.step (  );
+		L.step ( );
 	}
 	L.finalize( );
Index: /mpdm/CMakeLists.txt
===================================================================
--- /mpdm/CMakeLists.txt (revision 204)
+++ /mpdm/CMakeLists.txt (revision 213)
@@ -7,20 +7,21 @@
 ## Save all needed libraries in variable BdmLibs
 #SET(BdmLibs bdm itpp_debug)
-SET(BdmLibs bdm itpp)
-SET(ITppLibs itpp)
 
 IF(WIN32)
-	SET(BdmLibs ${BdmLibs} libacml_dll)
-	SET(ITppLibs ${ITppLibs} libacml_dll)
+	SET(AddLib libacml_dll)
 ENDIF(WIN32)
+
+# Define macro for testing a file
+MACRO(TEST FN)
+	add_executable (${FN} ${FN}.cpp)
+	target_link_libraries (${FN} debug itpp_debug)
+	target_link_libraries (${FN} optimized itpp)
+	target_link_libraries (${FN} bdm ${AddLib})
+ENDMACRO(TEST)
 
 # Add executable called "helloDemo" that is built from the source files
 # "demo.cxx" and "demo_b.cxx".  The extensions are automatically found.
-add_executable (merg_pred merg_pred.cpp)
-target_link_libraries (merg_pred ${BdmLibs})
-
-add_executable (merg_giw merg_giw.cpp)
-target_link_libraries (merg_giw ${BdmLibs})
-
-add_executable (merger_iter_cond merger_iter_cond.cpp)
-target_link_libraries (merger_iter_cond ${BdmLibs})
+TEST(merg_pred)
+TEST(merg_giw)
+TEST(merger_iter_cond)
+TEST(merg_2a)
Index: /mpdm/merg_2a.cpp
===================================================================
--- /mpdm/merg_2a.cpp (revision 211)
+++ /mpdm/merg_2a.cpp (revision 213)
@@ -10,30 +10,35 @@
 using std::endl;
 
+/*! The purpose of this experiment is to test merging fragmental pdfs between two participants sharing the same parameter, a. However, this parameter has a different role in each model.
+
+Lets assume that:
+ u -> |a,r| -> y -> |a,r| -> z
+*/
+
 int main() {
 	// Setup model
 	RV y ( "{y }" );
 	RV u ( "{u }" );
-	RV um = u; um.t(-1);
 	RV z ( "{z }" );
-	RV a ("{a }");
-	RV b ("{b }");
-	RV c ("{c }");
+	RV a ("{a }"); 
+	RV b ("{b }"); RV ab = a; ab.add(b);
+	RV c ("{c }"); RV ac = a; ac.add(c);
 	RV r ("{r }");
 	
 	double at = 1.5;
-	double bt = 0.8;
-	double ct = 0.50;
-	double sig = 0.10;
+	double bt = 0.5;
+	double ct = -0.5;
+	double rt = 0.30;
 	// Full system
 	vec thy =vec_2 ( at,bt ); //Simulated system - zero for constant term
-	vec thu =vec_2 ( at,ct ); //Simulated system - zero for constant term
-	vec Thy = concat ( thy, vec_1(sig) ); //Full parameter
-	vec Thu = concat ( thu, vec_1(sig) ); //Full parameter
+	vec thz =vec_2 ( at,ct ); //Simulated system - zero for constant term
+	vec Thy = concat ( thy, vec_1(rt) ); //Full parameter
+	vec Thz = concat ( thz, vec_1(rt) ); //Full parameter
 
 	//ARX constructor
 	mat V0 = 0.001*eye ( 3 ); V0 ( 0,0 ) = 1; //
 
-	ARX P1 ( concat ( a, concat(b,r) ), V0, -1 );
-	ARX P2 ( concat ( a, concat(c,r) ), V0, -1 );
+	ARX P1 ( concat ( ab,r ), V0, -1 );
+	ARX P2 ( concat ( ac,r ), V0, -1 );
 
 	//Test estimation
@@ -42,42 +47,40 @@
 
 	// Logging
-	dirfilelog L ( "exp/merg_2a",ndat );
-	int Li_Data = L.add ( RV ( "{Y U1 U2 }" ), "" );
-	int Li_LL   = L.add ( RV ( "{G M }" ), "LL" );
-	int Li_P1m   = L.add ( RV ( "{a b r }" ), "P1" );
-	int Li_P2m   = L.add ( RV ( "{a c r }" ), "P2" );
+	dirfilelog L ( "exp/merg_2a",3 );
+	int Li_Data = L.add ( RV ( "{U Y Z }" ), "" );
+//	int Li_LL   = L.add ( RV ( "{P1 P2 M1 M2 }" ), "LL" );
+	int Li_P1m   = L.add ( concat ( ab,r ), "P1m" );
+	int Li_P2m   = L.add ( concat ( ac,r ), "P2m" );
+	int Li_Mm   = L.add ( concat ( ab,concat(r,c) ), "Mm" );
 	L.init();
 
+	vec Ut ( ndat );
 	vec Yt ( ndat );
+	vec Zt ( ndat );
 	vec yt ( 1 );
 
 	//Proposal
-	enorm<ldmat> g0 ( a1 ); g0.set_parameters ( "1  ",mat ( "1" ) );
+	enorm<ldmat> g0 ( ab ); g0.set_parameters ( "1 1 ",mat ( "1 0; 0 1" ) );
 	egamma g1 ( r ); g1.set_parameters ( "2  ", "2" );
-	enorm<ldmat> g2 ( a2 ); g2.set_parameters ( "1  ",mat ( "1" ) );
+	enorm<ldmat> g2 ( c ); g2.set_parameters ( "1 ",mat ( "1" ) );
 
-	Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A ( 2 ) = &g2;
+	Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A(2) = &g2; 
 	eprod G0 ( A );
 
-	for ( t=0; t<ndat; t++ ) {
+	vec rgru(2);
+	vec rgry(2);
+	Yt(0) = 0.1;
+	Ut(0) = 0.0;
+	for ( t=1; t<ndat; t++ ) {
 		// True system
-		rgrg ( 0 ) = pow ( sin ( ( t/40.0 ) *pi ),3 );
-		rgrg ( 1 ) = pow ( cos ( ( t/40.0 +0.1 ) *pi ),3 );
-
-		Yt ( t ) = thg*rgrg + sqr * NorRNG();
+		Ut ( t ) = pow ( sin ( ( t/40.0 ) *pi ),3 );
+		rgru(0) = Ut(t); rgru(1) = Ut(t-1);
+		Yt ( t ) = thy*rgru + rt * NorRNG();
+		rgry(0) = Yt(t); rgry(1) = Yt(t-1);
+		Zt ( t ) = thz*rgry + rt * NorRNG();
 
 		// Bayes for all
-		P1.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 0 ) ) ) );
-		P2.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 1 ) ) ) );
-		PG.bayes ( concat ( Yt ( t ),rgrg ) );
-
-
-		// crippling PGk by substituting zeros.
-		/*	ldmat &Vk=const_cast<egiw*>(PGk._e())->_V();  //PG ldmat does not like 0!
-			mat fVk=PG._e()->_V().to_mat();
-			fVk(1,2) = 0.0;
-			fVk(2,1) = 0.0;
-			Vk = ldmat(fVk);
-		*/	//PGk is now krippled
+		P1.bayes ( concat ( Yt ( t ),rgru ) );
+		P2.bayes ( concat ( Zt ( t ),rgry ) );
 
 		// Merge estimates
@@ -89,31 +92,18 @@
 		//M2.set_parameters ( 100.0, 1000,3 ); //M2._Mix().set_method(QB);
 		M.merge ( &G0 );
-		//M2.merge ( &G0 );
 
 		//Likelihood
 		yt ( 0 ) = Yt ( t );
 
-//		LLs ( 0 ) = P1._e()->evalpdflog ( get_vec(Th, "1 2") );
-//		LLs ( 1 ) = P2._e()->evalpdflog ( get_vec(Th, "3 2")  );
-		LLs ( 0 ) = PG._e()->evalpdflog ( Th );
-		LLs ( 1 ) = M._Mix().logpred ( concat ( Thj, vec_1 ( 1.0 ) ) );
-//		LLs ( 3 ) = M2._Mix().logpred ( concat(Thj, vec_1(1.0)) );
-		L.logit ( Li_LL, LLs ); //log-normal
-
 		//Logger
-		L.logit ( Li_Data, vec_3 ( Yt ( t ), rgrg ( 0 ), rgrg ( 1 ) ) );
-		L.logit ( Li_P1m, P1._e()->mean() );
-		L.logit ( Li_P2m, P2._e()->mean() );
-		L.logit ( Li_Gm, PG._e()->mean() );
-		L.logit ( Li_Mm, M.mean() );
-
+		L.logit(Li_Data, vec_3(Ut(t), Yt(t), Zt(t)));
+		L.logit(Li_P1m, P1._e()->mean());
+		L.logit(Li_P2m, P2._e()->mean());
+		L.logit(Li_Mm, M.mean());
 		L.step ( );
 
-		cout << "Vg: " << PG._e()->_V().to_mat() <<endl;
-		vec mea = M.mean();
-		cout << "Ve: " << M.variance() <<endl;
 	}
 	L.finalize( );
-	L.itsave ( "merg_egiw.it" );
+	L.itsave ( "merg_2a.it" );
 	cout << endl;
 }
