/*!
  \file
  \brief TR 2525 file for testing Toy Problem of mpf for Covariance Estimation
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/



#include <estim/particles.h>
#include <estim/ekf_template.h>
#include <base/loggers.h>


#include "../pmsm.h"
#include "simulator.h"
#include "../sim_profiles.h"

using namespace bdm;

int main ( int argc, char* argv[] ) {
	const char *fname;
	if ( argc>1 ) {fname = argv[1];	}
	else { fname = "unitsteps.cfg"; }
	UIFile F ( fname );

	double h = 1e-6;
	int Nsimstep = 125;


	// Kalman filter
	int Ndat;
	int Npart;
	F.lookupValue ( "ndat", Ndat );
	F.lookupValue ( "Npart",Npart );
	shared_ptr<pdf> evolQ = UI::build<pdf>( F, "Qrw" );
	vec Qdiag;
	vec Rdiag;
	UI::get( Qdiag, F, "dQ" ); //( "1e-6 1e-6 0.001 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001
	UI::get( Rdiag, F, "dR" );// ( "1e-8 1e-8" ); //var(diff(xth)) = "0.034 0.034"

	// internal model

	shared_ptr<IMpmsm> fxu= new IMpmsm;
//                  Rs    Ls        dt       Fmag(Ypm)    kp   p    J     Bf(Mz)
	fxu->set_parameters ( 0.28, 0.003465, Nsimstep*h, 0.1989, 1.5 ,4.0, 0.04, 0.0 );
	// observation model
	shared_ptr<OMpmsm> hxu=new OMpmsm;

	vec mu0= "0.0 0.0 0.0 0.0";
	chmat Q ( Qdiag );
	chmat R ( Rdiag );
	EKFCh KFE ;
	KFE.set_parameters ( fxu,hxu,Q,R );
	KFE.set_statistics ( mu0, chmat ( zeros ( 4 ) ) );
	KFE.set_rv ( rx );
	KFE.validate();

	RV rQ ( "{Q }","4" );
	RV rU ("{u }","2");
	RV rY ("{y }","2");
	EKFCh_dQ KFEp ;
	KFEp.set_parameters ( fxu,hxu,Q,R );
	KFEp.set_statistics ( mu0, chmat ( zeros ( 4 ) ) );
	KFEp.set_rv(rx);
	KFEp.set_yrv(rY); 
	KFEp.set_rvc(concat(rU, rQ)); 
	KFEp.validate();

	MPF M;
	evolQ->set_rv(rQ);
	M.set_pf ( evolQ,Npart );
	M._pf().set_statistics(ones(Npart), euni(zeros(4),2*Qdiag));
	M.set_BM(KFEp);
	M.set_yrv ( rY   );
	M.set_rvc ( rU   );
	M.validate();

	shared_ptr<dirfilelog> L = UI::build<dirfilelog>( F, "logger" );// ( "exp/mpf_test",100 );
	int l_X = L->add_vector ( rx, "xt" );
	int l_D = L->add_vector ( concat ( ry,ru ), "" );
	int l_Q= L->add_vector ( rQ, "" );

	KFE.log_level[ logbounds] = true;
	KFE.log_register ( *L,"KF" );
	M.log_level[logbounds] = true;
	M.log_register ( *L,"M" );
	L->init();

	// SET SIMULATOR
	pmsmsim_set_parameters ( 0.28,0.003465,0.1989,0.0,4,1.5,0.04, 200., 3e-6, h );
	vec dt ( 2 );
	vec ut ( 2 );
	vec xt ( 4 );
	vec xtm=zeros ( 4 );
	double Ww=0.0;
	vec vecW;
	UI::get( vecW, F, "profile" );

	for ( int tK=1;tK<Ndat;tK++ ) {
		//Number of steps of a simulator for one step of Kalman
		for ( int ii=0; ii<Nsimstep;ii++ ) {
			//simulator
			sim_profile_vec01t ( Ww,vecW );
			pmsmsim_step ( Ww );
		};
		ut ( 0 ) = KalmanObs[4];
		ut ( 1 ) = KalmanObs[5];
		xt = fxu->eval ( xtm,ut ) + diag ( sqrt ( Qdiag ) ) *randn ( 4 );
		dt = hxu->eval ( xt,ut ) + diag(sqrt(Rdiag))*randn(2);
		xtm = xt;

		//Variances
		if ( tK==1000 )  Qdiag ( 0 ) *=10;
		if ( tK==2000 ) Qdiag ( 0 ) /=10;
		if ( tK==3000 )  Qdiag ( 1 ) *=10;
		if ( tK==4000 ) Qdiag ( 1 ) /=10;
		if ( tK==5000 )  Qdiag ( 2 ) *=10;
		if ( tK==6000 ) Qdiag ( 2 ) /=10;
		if ( tK==7000 )  Qdiag ( 3 ) *=10;
		if ( tK==8000 ) Qdiag ( 3 ) /=10;

		//estimator
		KFE.bayes (  dt,ut  );
		M.bayes ( dt,ut );

		L->log_vector ( l_X,xt );
		L->log_vector ( l_D,concat ( dt,ut ) );
		L->log_vector ( l_Q,Qdiag );

		KFE.log_write ( );
		M.log_write (  );
		L->step();
	}
	L->finalize();
	//Exit program:

	return 0;
}