#include <estim/arx.h>
#include <stat/libEF.h>
#include <stat/loggers.h>
#include <stat/merger.h>
using namespace itpp;

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

int main() {
	// Setup model
	RV y ( "{y }" );
	RV u1 ( "{u1 }" );
	RV u2 ( "{u2 }" );

	// Full system
	vec thg ( "0.7 1 1" ); //Simulated system
	//y=a y_t-1 + u1 + u2
	double sqr=0.1;
	int ord = 1;

	// Estimated systems ARX(2)
	RV thri ( "{thr_i }",vec_1 ( 2+1 ) );
	RV thrg ( "{thr_g }",vec_1 ( 3+1 ) );
	// Setup values

	//ARX constructor
	mat V0 = 0.001*eye ( thri.count() ); V0 ( 0,0 ) *= 10; //
	mat V0g = 0.001*eye ( thrg.count() ); V0g ( 0,0 ) *= 10; //
	double nu0 = ord+1;
	double frg = 0.99;

	ARX P1 ( thri, V0, nu0, frg );
	ARX P2 ( thri, V0, nu0, frg );
	ARX PG ( thrg, V0g, nu0, frg );
	//Test estimation
	int ndat = 10000;
	int t;

	// Logging
	dirfilelog L ( "exp/merg",ndat );

	int Eth1_log = L.add ( thri,"P1" );
	int Eth2_log = L.add ( thri,"P2" );
	int Ethg_log = L.add ( thrg,"PG" );
	int Data_log = L.add ( RV ( "{Y U1 U2 }" ), "" );
	int LL_log   = L.add ( RV ( "{1 2 G }" ), "LL" );

	L.init();

	vec Yt ( ndat );

	Yt.set_subvector ( 0,randn ( ord ) ); //initial values
	vec rgrg ( thg.length() );
	vec rgri ( thri.count() );

	for ( t=0; t<ndat; t++ ) {
		// True system
		if ( t>0 ) {
			rgrg ( 0 ) =Yt ( t-1 );
			rgrg ( 1 ) = pow(sin ( ( t/40.0 ) *pi ),3);
			rgrg ( 2 ) = pow(cos ( ( t/40.0 ) *pi ),3);

			Yt ( t ) = thg*rgrg + sqr * NorRNG();

			// 1st
			rgri ( 0 ) =Yt ( t );
			rgri ( 1 ) =Yt ( t-1 );
			rgri ( 2 ) =rgrg ( 1 );
			P1.bayes ( rgri );
			// 2nd
			rgri ( 2 ) =rgrg ( 2 );
			P2.bayes ( rgri );

			//Global
			PG.bayes ( concat ( Yt ( t ),rgrg ) );
			
			//Merger
		}
		L.logit ( Eth1_log, P1._epdf().mean() );
		L.logit ( Eth2_log, P2._epdf().mean() );
		L.logit ( Ethg_log, PG._epdf().mean() );
		L.logit ( Data_log, vec_3 ( Yt ( t ), rgrg ( 1 ), rgrg ( 2 ) ) );
		L.logit ( LL_log, vec_3 ( P1._ll(), P2._ll(), PG._ll() ) );
		L.step (  );
	}
	L.finalize( );
	L.itsave ( "merg.it" );
}