/*
  \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"
#include "simulator.h"

using namespace itpp;
//!Extended Kalman filter with unknown \c Q
class EKF_unQ : public EKFCh , public BMcond {
public:
	//! Default constructor
	EKF_unQ ( RV rx, RV ry,RV ru,RV rQ ) :EKFCh ( rx,ry,ru ),BMcond ( rQ ) {};
	void condition ( const vec &Q0 ) {
		Q.setD ( Q0,0 );
		//from EKF
		preA.set_submatrix ( dimy+dimx,dimy,Q._Ch() );
	}; 
};


int main() {
	// Kalman filter
	int Ndat = 10000;
	double h = 1e-6;
	int Nsimstep = 125;
	int Npart = 100;
	
	// 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 ( "0.05 0.05 0.001 0.001" ); //zdenek: 0.05 0.05 0.001 0.001
	vec Rdiag ( "0.03 0.03" ); //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 ( 1000*ones ( 4 ) ) );

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

	mgamma evolQ ( rQ,rQ );
	MPF<EKF_unQ> M ( rx,rQ,evolQ,evolQ,Npart,KFEp );
	// initialize
	evolQ.set_parameters ( 100.0 ); //sigma = 1/10 mu
	evolQ.condition ( "0.05 0.05 0.001 0.001" ); //Zdenek default
	epdf& pfinit=evolQ._epdf();
	M.set_est ( pfinit );
	evolQ.set_parameters ( 1000.0 ); 

	//

	epdf& KFEep = KFE._epdf();
	epdf& Mep = M._epdf();

	mat Xt=zeros ( 9,Ndat ); //true state from simulator
	mat Dt=zeros ( 4,Ndat ); //observation
	mat XtE=zeros ( 4,Ndat );
	mat XtM=zeros ( 8,Ndat ); //Q + x

	// SET SIMULATOR
	pmsmsim_set_parameters ( 0.28,0.003465,0.1989,0.0,4,1.5,0.04, 200., 3e-6, h );
	double Ww=0.0;
	static int k_rampa=1;
	static long k_rampa_tmp=0;
	vec dt ( 2 );
	vec dtVS =zeros( 2 );
	vec xtVS =zeros(4);
	vec et ( 4 );
	vec wt ( 2 );
	vec ut ( 2 );
	mat XtV=zeros ( 4,Ndat );

	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
			Ww+=k_rampa*2.*M_PI*2e-4;    //1000Hz/s
			if ( Ww>2.*M_PI*150. ) {
				Ww=2.*M_PI*150.;
				if ( k_rampa_tmp<500000 ) k_rampa_tmp++;
				else {k_rampa=-1;k_rampa_tmp=0;}
			};
			if ( Ww<-2.*M_PI*150. ) Ww=-2.*M_PI*150.; /* */

			pmsmsim_step ( Ww );
		};
		// collect data
		ut ( 0 ) = KalmanObs[0];
		ut ( 1 ) = KalmanObs[1];
		dt ( 0 ) = KalmanObs[2];
		dt ( 1 ) = KalmanObs[3];
		
		// My own simulator for testing : Asuming ut is OK
		NorRNG.sample_vector ( 4,et );
		NorRNG.sample_vector ( 2,wt );
		xtVS = fxu.eval ( xtVS,ut ) + Q.sqrt_mult ( et );
		dtVS = hxu.eval ( xtVS,ut ) + R.sqrt_mult ( wt );

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

		Xt.set_col ( tK,vec ( x,9 ) ); //vec from C-array
		Dt.set_col ( tK, concat ( dt,ut ) );
		XtE.set_col ( tK,KFEep.mean() );
		XtM.set_col ( tK,Mep.mean() );
		XtV.set_col ( tK,xtVS );
	}

	it_file fou ( "pmsm_sim.it" );

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

}