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

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

#include "iopom.h"

using namespace itpp;
//!Extended Kalman filter with unknown \c Q

void set_simulator_t ( double &Ww ) {

	if ( t>0.0002 ) x[8]=1.2;  // 1A //0.2ZP
	if ( t>0.4 ) x[8]=10.8; // 9A
	if ( t>0.6 ) x[8]=25.2;  // 21A

	if ( t>0.7 ) Ww=2.*M_PI*10.;
	if ( t>1.0 ) x[8]=1.2;  // 1A
	if ( t>1.2 ) x[8]=10.8; // 9A
	if ( t>1.4 ) x[8]=25.2;  // 21A

	if ( t>1.6 ) Ww=2.*M_PI*50.;
	if ( t>1.9 ) x[8]=1.2;  // 1A
	if ( t>2.1 ) x[8]=10.8; // 9A
	if ( t>2.3 ) x[8]=25.2;  // 21A

	if ( t>2.5 ) Ww=2.*M_PI*100;
	if ( t>2.8 ) x[8]=1.2;  // 1A
	if ( t>3.0 ) x[8]=10.8; // 9A
	if ( t>3.2 ) x[8]=25.2;  // 21A

	if ( t>3.4 ) Ww=2.*M_PI*150;
	if ( t>3.7 ) x[8]=1.2;  // 1A
	if ( t>3.9 ) x[8]=10.8; // 9A
	if ( t>4.1 ) x[8]=25.2;  // 21A

	if ( t>4.3 ) Ww=2.*M_PI*0;
	if ( t>4.8 ) x[8]=-1.2;  // 1A
	if ( t>5.0 ) x[8]=-10.8; // 9A
	if ( t>5.2 ) x[8]=-25.2;  // 21A

	if ( t>5.4 ) Ww=2.*M_PI* ( -10. );
	if ( t>5.7 ) x[8]=-1.2;  // 1A
	if ( t>5.9 ) x[8]=-10.8; // 9A
	if ( t>6.1 ) x[8]=-25.2;  // 21A

	if ( t>6.3 ) Ww=2.*M_PI* ( -50. );
	if ( t>6.7 ) x[8]=-1.2;  // 1A
	if ( t>6.9 ) x[8]=-10.8; // 9A
	if ( t>7.1 ) x[8]=-25.2;  // 21A

	if ( t>7.3 ) Ww=2.*M_PI* ( -100. );
	if ( t>7.7 ) x[8]=-1.2;  // 1A
	if ( t>7.9 ) x[8]=-10.8; // 9A
	if ( t>8.1 ) x[8]=-25.2;  // 21A
	if ( t>8.3 ) x[8]=10.8; // 9A
	if ( t>8.5 ) x[8]=25.2;  // 21A

        x[8]=0.0;
}

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

	mat Xt=zeros ( Ndat ,9 ); //true state from simulator
	mat XtM=zeros ( Ndat ,4 ); //true state from simulator
	mat XtF=zeros ( Ndat ,4 ); //true state from simulator
	mat XtO=zeros ( Ndat ,4 ); //true state from simulator
	mat Dt=zeros ( Ndat,4 ); //observation
	mat Qrec=zeros ( Ndat, 4*4 ); //full covariance matrix

	// 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;
	vec dt ( 2 );
	vec ut ( 2 );
	vec xtm=zeros ( 4 );
	vec xdif=zeros ( 4 );
	vec xt ( 4 );
	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 );
	OMpmsm hxu;
	mat Qt=zeros ( 4,4 );

	// ESTIMATORS
	vec mu0= "0.0 0.0 0.0 0.0";
	vec Qdiag ( "62 66 454 0.03" ); //zdenek: 0.01 0.01 0.0001 0.0001
	vec Rdiag ( "100 100" ); //var(diff(xth)) = "0.034 0.034"
	mat Q =diag( Qdiag );
	mat R =diag ( Rdiag );
	EKFfull Efix ( rx,ry,ru );
	Efix.set_est ( mu0, 1*eye ( 4 )  );
	Efix.set_parameters ( &fxu,&hxu,Q,R);

	EKFfull Eop ( rx,ry,ru );
	Eop.set_est ( mu0, 1*eye ( 4 ) );
	Eop.set_parameters ( &fxu,&hxu,Q,R);

	epdf& Efix_ep = Efix._epdf();
	epdf& Eop_ep = Eop._epdf();


	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++ ) {
			set_simulator_t ( Ww );
			pmsmsim_step ( Ww );
		};
		// simulation via deterministic model
		ut ( 0 ) = KalmanObs[0];
		ut ( 1 ) = KalmanObs[1];
		dt ( 0 ) = KalmanObs[2];
		dt ( 1 ) = KalmanObs[3];

		xt = fxu.eval ( xtm,ut );

		//Results:
		// in xt we have simulaton according to the model
		// in x we have "reality"

		xtm ( 0 ) =x[0];xtm ( 1 ) =x[1];xtm ( 2 ) =x[2];xtm ( 3 ) =x[3];
		xdif=xtm-xt;
		if ( xdif ( 3 ) >pi ) xdif ( 3 )-=2*pi;
		if ( xdif ( 3 ) <-pi ) xdif ( 3 ) +=2*pi;

		//Rt = 0.9*Rt + xdif^2
		Qt*=0.0;
		Qt += outer_product ( xdif,xdif ); //(x-xt)^2
		Qrec.set_row ( tK, vec ( Qt._data(),16 ) );

		//ESTIMATE
		Efix.bayes(concat(dt,ut));
		//
		Eop.set_parameters ( &fxu,&hxu,(Qt+1e-16*eye(4)),(1*eye(2)));
		Eop.bayes(concat(dt,ut));
		
		Xt.set_row ( tK,vec ( x,9 ) ); //vec from C-array
		XtM.set_row ( tK,xt ); //vec from C-array
		XtF.set_row ( tK,Efix_ep.mean() ); 
		XtO.set_row ( tK,Eop_ep.mean() ); 
		Dt.set_row ( tK, concat ( dt,ut ) );
	}

	char tmpstr[200];
	sprintf ( tmpstr,"%s/%s","variance/","format" );
	ofstream  form ( tmpstr );
	form << "# Experimetal header file"<< endl;
	dirfile_write ( form,"variance/",Xt,"X","{isa isb om th }" );
	dirfile_write ( form,"variance/",XtM,"Xsim","{isa isb om th }" );
	dirfile_write ( form,"variance/",XtO,"XO","{isa isb om th }" );
	dirfile_write ( form,"variance/",XtF,"XF","{isa isb om th }" );
	dirfile_write ( form,"variance/",Dt,"D","{isa isb usa usb }" );
	dirfile_write ( form,"variance",Qrec,"Q","{ }" );

	return 0;
}