/*
  \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"

#include "../simulator_zdenek/ekf_example/ekf_obj.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;
	int Npart = 100;
	

	vec mu0= "0.0 0.0 0.0 0.0";
	vec Qdiag ( "0.05 0.05 0.002 0.001" ); //zdenek: 0.01 0.01 0.0001 0.0001
	vec Rdiag ( "0.05 0.05" ); //var(diff(xth)) = "0.034 0.034"
	chmat Q ( Qdiag );
	chmat R ( Rdiag );
	
	RV rQ ( "100","{Q}","4","0" );
	EKFfixed KFE ( rx, rQ);
	KFE.init_ekf ( Nsimstep*h);


///////////// Particles
	mgamma_fix evolQ ( rQ,rQ );
	MPF<EKFfixed> M ( rx,rQ,evolQ,evolQ,Npart,KFE );
	// initialize
	evolQ.set_parameters ( 10.0 ,Qdiag, 1.0); //sigma = 1/10 mu
	evolQ.condition ( Qdiag ); //Zdenek default
	epdf& pfinit=evolQ._epdf();
	M.set_est ( pfinit );
	evolQ.set_parameters ( 500000.0, Qdiag, 0.9999 );


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


	mat Xt=zeros ( Ndat ,9 ); //true state from simulator
	mat Dt=zeros ( Ndat,4 ); //observation
	mat XtE=zeros ( Ndat, 4 );
	mat XtM=zeros ( Ndat,4+4 ); //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 ut ( 2 );

	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.; /* */
//			set_simulator_t(Ww);
			pmsmsim_step ( Ww );
		};
		// collect data
		ut ( 0 ) = KalmanObs[0];
		ut ( 1 ) = KalmanObs[1];
		dt ( 0 ) = KalmanObs[2];
		dt ( 1 ) = KalmanObs[3];
		//estimator
		KFE.bayes ( concat ( dt,ut ) );
		M.bayes ( concat ( dt,ut ) );
		
		Xt.set_row ( tK,vec ( x,9 ) ); //vec from C-array
		Dt.set_row ( tK, concat ( dt,ut));
		XtE.set_row ( tK,KFEep.mean() );
		XtM.set_row ( tK,Mep.mean() );
	}

	char tmpstr[200];
	sprintf(tmpstr,"%s/%s","herez/","format");
	ofstream  form(tmpstr);
	form << "# Experimetal header file"<< endl;
	dirfile_write(form,"herez/",Xt,"X","{isa isb om th }" );
	dirfile_write ( form,"herez",XtM,"XtM","{q1 q2 q3 q4 isa isb om th }" );
	dirfile_write ( form,"herez",XtE,"XE","{isa isb om th }" );
	dirfile_write ( form,"herez",Dt,"Dt","{isa isb ua ub }" );

	return 0;
}