#include <itpp/itbase.h>
#include <stat/libEF.h>
#include <estim/merger.h>

using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

int main() {

	RNG_randomize();

	RV y ( "{y }","1" );
	RV u1 ( "{u1 }","1" );
	RV u2 ( "{u2 }","1" );

	RV all = y;
	all.add ( u1 );
	all.add ( u2 );

	mlnorm<fsqmat> f1 ( y , u1 );
	mlnorm<fsqmat> f2 ( y , u2 );

	//Differneces in constant term are essential
	f1.set_parameters ( "0.4","1",mat ( "0.04" ) );
	f2.set_parameters ( "0.2","-1",mat ( "0.08" ) );

	Array<mpdf* > A ( 2 );
	A ( 0 ) =&f1;
	A ( 1 ) =&f2;

	merger M ( A );
	enorm<fsqmat> g0 ( all );
	g0.set_parameters ( vec ( "1 -1 1" ),3*eye ( 3 ));// +1*ones ( 3,3 ) );

	M.set_parameters ( 1.2,1000,1 );

	int Ntrials=100;
	vec A1s ( Ntrials );
	vec A2s ( Ntrials );
	vec R1s ( Ntrials );
	vec R2s ( Ntrials );
	vec C1s ( Ntrials );
	vec C2s ( Ntrials );

	for ( int it=0;it<Ntrials;it++ ) {
		M.merge ( &g0 );

		MixEF &MM = M._Mix();
		epdf* MP = MM._Coms ( 0 )->predictor ( all );

		RV yu1 = y; yu1.add ( u1 );
		RV yu2 = y; yu2.add ( u2 );
		enorm<ldmat>* P1m= ( enorm<ldmat>* ) MP->marginal ( yu1 );
		enorm<ldmat>* P2m= ( enorm<ldmat>* ) MP->marginal ( yu2 );
		mlnorm<ldmat>* P1c= ( mlnorm<ldmat>* ) ( P1m->condition ( y ) );
		mlnorm<ldmat>* P2c= ( mlnorm<ldmat>* ) ( P2m->condition ( y ) );

		A1s(it) = P1c->_A()(0,0);
		A2s(it) = P2c->_A()(0,0);
		R1s(it) = P1c->_R()(0,0);
		R2s(it) = P2c->_R()(0,0);
		C1s(it) = P1c->_mu_const()(0);
		C2s(it) = P2c->_mu_const()(0);

	}
	double A1mean = sum(A1s)/Ntrials;
	double A2mean = sum(A2s)/Ntrials;
	double C1mean = sum(C1s)/Ntrials;
	double C2mean = sum(C2s)/Ntrials;
	double R1mean = sum(R1s)/Ntrials;
	double R2mean = sum(R2s)/Ntrials;
	cout << "A1: " << A1mean << " +- " << 2*(sum_sqr(A1s)/Ntrials-A1mean*A1mean) <<endl;
	cout << "A2: " << A2mean << " +- " << 2*(sum_sqr(A2s)/Ntrials-A2mean*A2mean) <<endl;
	cout << "C1: " << C1mean << " +- " << 2*(sum_sqr(C1s)/Ntrials-C1mean*C1mean) <<endl;
	cout << "C2: " << C2mean << " +- " << 2*(sum_sqr(C2s)/Ntrials-C2mean*C2mean) <<endl;
	cout << "R1: " << R1mean << " +- " << 2*(sum_sqr(R1s)/Ntrials-R1mean*R1mean) <<endl;
	cout << "R2: " << R2mean << " +- " << 2*(sum_sqr(R2s)/Ntrials-R2mean*R2mean) <<endl;
}