#include <itpp/itbase.h>
#include <estim/libKF.h>

using namespace itpp;

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

int main() {


	// Klaman filter
	mat A, B,C,D,R,Q,P0;
	vec mu0;
	mat Mu0;;
	// input from Matlab
	it_file fin( "testKF.it" );

	mat Dt, Xt,Xt2,XtE,Xtf;
	int Ndat;

	bool xxx= fin.seek( "d" );
	if (!xxx){ it_error("testKF.it not found");}
	fin >>Dt;
	fin.seek( "A" ); 
	fin >> A;
	fin.seek( "B" ); 
	fin >> B;
	fin.seek( "C" ); 
	fin >> C;
	fin.seek( "D" ); 
	fin >> D;
	fin.seek( "R" ); 
	fin >> R;
	fin.seek( "Q" ); fin >> Q;
	fin.seek( "P0" ); fin >> P0;
	fin.seek( "mu0" ); fin >> Mu0; 
	mu0=Mu0.get_col(0);
	
	Ndat = Dt.cols();
	Xt=zeros( 2,Ndat );
	Xt2=zeros( 2,Ndat );
	Xtf=zeros( 2,Ndat );
	XtE=zeros( 2,Ndat );

//	cout << KF;
	RV rx("1","{x}","2","0");
	RV ru("2","{u}","1","0");
	RV ry("3","{y}","1","0");
	//
	Kalman<ldmat> KF(rx,ry,ru);
	KF.set_parameters(A,B,C,D,ldmat(R),ldmat(Q));
	KF.set_est(mu0,ldmat(P0) );
	//
	Kalman<fsqmat> KFf(rx,ry,ru);
	KFf.set_parameters(A,B,C,D,fsqmat(R),fsqmat(Q));
	KFf.set_est(mu0,fsqmat(P0) );
	//
	KalmanFull KF2( A,B,C,D,R,Q,P0,mu0 );
	//
	bilinfn fxu(rx,ru,A,B);
	bilinfn hxu(rx,ru,C,D);
	EKF<ldmat> KFE(rx,ry,ru);
	KFE.set_parameters(&fxu,&hxu,Q,R);
	KFE.set_est(mu0,P0);

	epdf& KFep = KF._epdf();
	epdf& KFfep = KFf._epdf();
	epdf& KFEep = KFE._epdf();
	
	Xt.set_col( 0,KFep.mean() );
	Xtf.set_col( 0,KFfep.mean() );
	Xt2.set_col( 0,KF2.mu );
	XtE.set_col( 0,KFEep.mean() );
	for ( int t=1;t<Ndat;t++ ) {
		KFf.bayes( Dt.get_col( t ));
		KF.bayes( Dt.get_col( t ));
		KF2.bayes( Dt.get_col( t ));
		KFE.bayes( Dt.get_col( t ));
		Xt.set_col( t,KFep.mean() );
		Xtf.set_col( t,KFfep.mean() );
		Xt2.set_col(t,KF2.mu);
		XtE.set_col( t,KFEep.mean() );
	}

	it_file fou( "testKF_res.it" );
	fou << Name("xth") << Xt;
	fou << Name("xthf") << Xtf;
	fou << Name("xth2") << Xt2;
	fou << Name("xthE") << XtE;
	//Exit program:
	return 0;

}