/*
  \file
  \brief Models for synchronous electric drive using IT++ and BDM
  \author Vaclav Smidl.

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

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

#include <itpp/itbase.h>
#include <estim/libKF.h>
#include <estim/libPF.h>
#include <stat/libFN.h>

#include "pmsm.h"

using namespace itpp;

//!Extended Kalman filter with unknown \c Q
class EKF_unQ : public EKF<ldmat> , public BMcond {
public:
	//! Default constructor
	EKF_unQ (RV rx, RV ry,RV ru,RV rQ ) :EKF<ldmat> ( rx,ry,ru ),BMcond ( rQ ) {};
	void condition ( const vec &Q0 ) {Q.setD ( Q0,0 ); }; //set first two diag(Q)
};

int main() {
	// Kalman filter
	int Ndat = 10000;

//	cout << KF;
	// internal model
	IMpmsm fxu;
	//                  Rs    Ls        dt       Fmag(Ypm) kp   p    J     Bf(Mz)
	fxu.set_parameters ( 0.28, 0.003465, 20*1e-6, 0.1989,   1.5 ,4.0, 0.04, 0.0 );
	// observation model
	OMpmsm hxu;

	vec mu0= "100 100 100 1";
	vec Qdiag ( "0.1 0.1 0.01 0.00001" );
	vec Rdiag ( "0.02 0.02" );
	vec vQ = "0.01:0.01:100";

	ldmat Q = ldmat ( Qdiag );
	ldmat R = ldmat ( Rdiag );

	RV rQ ( "100","{Q}","2","0" );
	EKF_unQ KFE ( rx,ry,ru,rQ );
	KFE.set_parameters ( &fxu,&hxu,Q,R );
	KFE.set_est ( mu0, ldmat ( 1000*ones ( 4 ) ) );

	mgamma evolQ ( rQ,rQ );
	evolQ.set_parameters ( 10000.0 ); //sigma = 1/10 mu

	MPF<EKF_unQ> M ( rx,rQ,evolQ,evolQ,100,KFE );

	epdf& KFEep = KFE._epdf();
	epdf& Mep = M._epdf();
	// initialize
	evolQ.condition ( "0.5 0.5" );
	epdf& pfinit=evolQ._epdf();
	M.set_est ( pfinit );

	//simulator values
	vec dt ( 2 );
	vec wt ( 2 );
	vec ut ( 2 );
	vec xt=mu0;
	vec et ( 4 );

	mat Xt=zeros ( 4,Ndat );
	mat XtE=zeros ( 4,Ndat );
	mat XtM=zeros ( 6,Ndat );
	Xt.set_col ( 0,KFEep.mean() );
	XtM.set_col ( 0,Mep.mean() );

	for ( int t=1;t<Ndat;t++ ) {
		//simulator
		UniRNG.sample_vector ( 2,wt );

		if ( rem ( t,500 ) <200 ) ut = rem ( t,500 ) *ones ( 2 );
		else
			ut=zeros ( 2 );

		NorRNG.sample_vector ( 4,et );
		NorRNG.sample_vector ( 2,wt );
		xt = fxu.eval ( xt,ut ) + Q.sqrt_mult ( et );
		dt = hxu.eval ( xt,ut ) + R.sqrt_mult ( wt );

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

		Xt.set_col ( t,xt );
		XtE.set_col ( t,KFEep.mean() );
		XtM.set_col ( t,Mep.mean() );
	}

	it_file fou ( "pmsm.it" );

	fou << Name ( "xth" ) << Xt;
	fou << Name ( "xthE" ) << XtE;
	fou << Name ( "xthM" ) << XtM;
	//Exit program:
	return 0;

}