/*!
  \file
  \brief Voltage U is multiplied by an unknown weight w which is estimated by MPF
  \author Vaclav Smidl.

  \ingroup PMSM

  -----------------------------------
  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 "sim_profiles.h"

using namespace itpp;

//!Extended Kalman filter with unknown \c Q
class EKFCh_cond : public EKFCh , public BMcond {
public:
	//! Default constructor
	EKFCh_cond ( RV rx, RV ry,RV ru,RV rC ) :EKFCh ( rx,ry,ru ),BMcond ( rC ) {};
	void condition ( const vec &val ) {
		pfxu->condition ( val );
	};
};

class IMpmsm_w :  public IMpmsm {
protected:
	double w;
public:
	IMpmsm_w() :IMpmsm(),w ( 1.0 ) {};
	//! Set mechanical and electrical variables

	void condition ( const vec &val ) {w = val ( 0 );}
	vec eval ( const vec &x0, const vec &u0 ) {
		// last state
		double iam = x0 ( 0 );
		double ibm = x0 ( 1 );
		double omm = x0 ( 2 );
		double thm = x0 ( 3 );
		double uam = u0 ( 0 );
		double ubm = u0 ( 1 );

		vec xk=zeros ( 4 );
		//ia
		xk ( 0 ) = ( 1.0- Rs/Ls*dt ) * iam + Ypm/Ls*dt*omm * sin ( thm ) + w* uam*dt/Ls;
		//ib
		xk ( 1 ) = ( 1.0- Rs/Ls*dt ) * ibm - Ypm/Ls*dt*omm * cos ( thm ) + w* ubm*dt/Ls;
		//om
		xk ( 2 ) = omm + kp*p*p * Ypm/J*dt* ( ibm * cos ( thm )-iam * sin ( thm ) ) - p/J*dt*Mz;
		//th
		xk ( 3 ) = thm + omm*dt; // <0..2pi>
		if ( xk ( 3 ) >pi ) xk ( 3 )-=2*pi;
		if ( xk ( 3 ) <-pi ) xk ( 3 ) +=2*pi;
		return xk;
	}

};


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

	mat Rnoise = randn ( 2,Ndat );

	// internal model
	IMpmsm fxu0;
	IMpmsm_w 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 );
	fxu0.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.1 0.1 0.001 0.000001" ); //zdenek: 0.01 0.01 0.0001 0.0001
	vec Rdiag ( "0.1 0.1" ); //var(diff(xth)) = "0.034 0.034"
	chmat Q ( Qdiag );
	chmat R ( Rdiag );
	EKFCh KFE ( rx,ry,ru );
	KFE.set_parameters ( &fxu0,&hxu,Q,R );
	KFE.set_est ( mu0, chmat ( ones ( 4 ) ) );

	RV rW ( "{w }" );
	EKFCh_cond KFEp ( rx,ry,ru,rW );
	KFEp.set_parameters ( &fxu,&hxu,Q,R );
	KFEp.set_est ( mu0, chmat ( ones ( 4 ) ) );

	mgamma_fix evolW ( rW,rW );
	MPF<EKFCh_cond> M ( rx,rW,evolW,evolW,Npart,KFEp );
	// initialize
	vec W0="0.5";
	evolW.set_parameters ( 10.0, W0, 1.0 ); //sigma = 1/10 mu
	evolW.condition ( W0 );
	epdf& pfinit=evolW._epdf();
	M.set_est ( pfinit );
	evolW.set_parameters ( 100.0, W0, 0.99 ); //sigma = 1/10 mu

	mat Xt=zeros ( Ndat ,4 ); //true state from simulator
	mat Dt=zeros ( Ndat,2+2 ); //observation
	mat XtE=zeros ( Ndat, 4 );
	mat Qtr=zeros ( Ndat, 4 );
	mat XtM=zeros ( Ndat,1+4 ); //W + x
	mat XtMTh=zeros ( Ndat,1 ); //W + x

	// 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 8 4 2 0 -4 -9 -16 -4 0 0 0";
	vecW*=10.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[0];
		ut ( 1 ) = KalmanObs[1];
		dt ( 0 ) = KalmanObs[2]+0.3*Rnoise ( 0,t );
		dt ( 1 ) = KalmanObs[3]+0.3*Rnoise ( 1,t );
		xt = vec ( x,4 );

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

		Xt.set_row ( tK, xt ); //vec from C-array
		Dt.set_row ( tK, concat ( dt,ut ) );
		Qtr.set_row ( tK, Qdiag );
		XtE.set_row ( tK,KFE._e()->mean() );
		XtM.set_row ( tK,M._e()->mean() );
		{
			double sumSin=0.0;
			double sumCos=0.0;
			vec mea ( 4 );
			vec* _w;

			for ( int p=0; p<Npart;p++ ) {
				mea = M._BM ( p )->_e()->mean();
				_w = M.__w();
				sumSin += ( *_w ) ( p ) *sin ( mea ( 3 ) );
				sumCos += ( *_w ) ( p ) *cos ( mea ( 3 ) );
			}
			double Th = asin ( sumSin );
			if ( sumCos<0 ) {
				if ( sumSin>0 ) {
					Th = M_PI-Th;
				}
				else {
					Th = -M_PI-Th;
				}
			}

			XtMTh.set_row ( tK,vec_1 ( Th ) );
		}
	}

	it_file fou ( "mpf_u_weight.it" );

	fou << Name ( "xth" ) << Xt;
	fou << Name ( "Dt" ) << Dt;
	fou << Name ( "Qtr" ) << Qtr;
	fou << Name ( "xthE" ) << XtE;
	fou << Name ( "xthM" ) << XtM;
	fou << Name ( "xthMTh" ) << XtMTh;
	//Exit program:

	return 0;
}