#include <estim/arx.h>
#include <stat/libEF.h>
#include <stat/loggers.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("1","{y }","1","0");
	RV u1("2","{u1 }","1","0");
	RV u2("3","{u2 }","1","0");
	
	// 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("4","{thr }",vec_1(2+1),"0");
	RV thrg("5","{thr }",vec_1(3+1),"0");
	// 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.94;
	
	ARX P1(thri, V0, nu0, frg);
	ARX P2(thri, V0, nu0, frg);
	ARX PG(thrg, V0g, nu0, frg);
	//Test estimation
	int ndat = 1000;
	int t;
	
	// Logging
	dirfilelog L("exp/merg",1000);

	int Eth1_log = L.add(thri,"P1");
	int Eth2_log = L.add(thri,"P2");
	int Data_log = L.add(RV("10","{Y U1 U2 }","3","0"), "");
	int LL_log   = L.add(RV("11","{1 2 G }","3","0"), "LL");
	
	vec Yt(ndat);
	
	Yt.set_subvector(0,randn(ord)); //initial values
	vec rgrg(thg.length());
	vec rgri(thri.count());
	
	for (t=2; t<ndat; t++) {
		// True system
		rgrg(0)=Yt(t-1); 
		rgrg(1)= sin(t/(40)*pi);
		rgrg(2)= cos(t/(40)*pi);
		
		Yt(t) = thg*rgrg + sqr * NorRNG();
		
		// 1st
		rgri(0)=Yt(t);
		rgri(1)=rgrg(1);
		P1.bayes(rgri);
		// 2nd
		rgri(1)=rgrg(2);
		P2.bayes(rgri);
		
		//Global
		PG.bayes(concat(Yt(t),rgrg));
		L.logit(Eth1_log, P1._epdf().mean());
		L.logit(Eth2_log, P2._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(false);
	}
	L.step(true);
	
}