/*!
  \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/libKF.h>
#include <estim/libPF.h>
#include <estim/ekf_templ.h>
#include <stat/libFN.h>

#include <stat/loggers.h>

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

using namespace bdm;

int main() {
	// Kalman filter
	int Ndat = 9000;
	double h = 1e-6;
	int Nsimstep = 125;
	int Npart = 200;

	// internal model
	IMpmsm fxu;
	//                  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
	OMpmsm hxu;

	vec mu0= "0.0 0.0 0.0 0.0";
	vec Qdiag ( "1e-6 1e-6 0.001 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001
	vec Rdiag ( "1e-8 1e-8" ); //var(diff(xth)) = "0.034 0.034"
	chmat Q ( Qdiag );
	chmat R ( Rdiag );
	EKFCh KFE ( rx,ry,ru );
	KFE.set_parameters ( &fxu,&hxu,Q,R );
	KFE.set_est ( mu0, chmat ( zeros ( 4 ) ) );

	RV rQ ( "{Q }","4" );
	EKFCh_unQ KFEp ( rx,ry,ru,rQ );
	KFEp.set_parameters ( &fxu,&hxu,Q,R );
	KFEp.set_est ( mu0, chmat ( zeros ( 4 ) ) );

	//mgamma_fix evolQ ( rQ,rQ );
	migamma_fix evolQ ( rQ,rQ );
	MPF<EKFCh_unQ> M ( rx,rQ,evolQ,evolQ,Npart,KFEp );
	// initialize
	evolQ.set_parameters ( 0.1, 10*Qdiag, 1.0); //sigma = 1/10 mu
	evolQ.condition (10*Qdiag ); //Zdenek default
	M.set_est ( *evolQ._e() );
	evolQ.set_parameters ( 0.10, 10*Qdiag,0.9999); //sigma = 1/10 mu
	//

	const epdf& KFEep = KFE.posterior();
	const epdf& Mep = M.posterior();

	dirfilelog L("exp/mpf_test",100);
	int l_X = L.add(rx, "xt");
	int l_D = L.add(concat(ry,ru), "");
	int l_XE= L.add(rx, "xtE");
	int l_XM= L.add(concat(rQ,rx), "xtM");
	int l_VE= L.add(rx, "VE");
	int l_VM= L.add(concat(rQ,rx), "VM");
	int l_Q= L.add(rQ, "");
	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="1 2 4 9 4 2 0 -4 -9 -16 -4 0 0";

	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);
		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 ( concat ( dt,ut ) );
		M.bayes ( concat ( dt,ut ) );

		L.logit(l_X,xt);
		L.logit(l_D,concat(dt,ut));
		L.logit(l_XE,KFEep.mean());
		L.logit(l_XM,Mep.mean());
		L.logit(l_VE,KFEep.variance());
		L.logit(l_VM,Mep.variance());
		L.logit(l_Q,Qdiag);
		L.step();
	}
	L.finalize();
	//Exit program:

	return 0;
}