/*!
  \file
  \brief Bayesian Filtering for linear Gaussian models (Kalman Filter) and extensions
  \author Vaclav Smidl.

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

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

#ifndef EKFfix_H
#define EKFfix_H


#include <estim/kalman.h>
#include "fixed.h"
#include "matrix.h"
#include "matrix_vs.h"
#include "reference_Q15.h"
#include "parametry_motoru.h"
#include "mpf_double.h"

using namespace bdm;

double minQ(double Q);

void mat_to_int16(const imat &M, int16 *I);
void vec_to_int16(const ivec &v, int16 *I);
void UDtof(const mat &U, const vec &D, imat &Uf, ivec &Df, const vec &xref);

#ifdef XXX
/*!
\brief Extended Kalman Filter with full matrices in fixed point16 arithmetic

An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean.
*/
class EKFfixed : public BM {
public:
void init_ekf(double Tv);
void ekf(double ux, double uy, double isxd, double isyd);

/* Declaration of local functions */
void prediction(int16 *ux);
void correction(void);
void update_psi(void);

/* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/
 int16 Q[16]; /* matrix [4,4] */
 int16 R[4]; /* matrix [2,2] */

 int16 x_est[4];
 int16 x_pred[4];
 int16 P_pred[16]; /* matrix [4,4] */
 int16 P_est[16]; /* matrix [4,4] */
 int16 Y_mes[2];
 int16 ukalm[2];
 int16 Kalm[8]; /* matrix [5,2] */

 int16 PSI[16]; /* matrix [4,4] */

 int16 temp15a[16];

 int16 cA, cB, cC, cG, cH;  // cD, cE, cF, cI ... nepouzivane

 int32 temp30a[4]; /* matrix [2,2] - temporary matrix for inversion */
 enorm<fsqmat> E;
 mat Ry;
  
public:
	//! Default constructor
	EKFfixed ():BM(),E(),Ry(2,2){
		int16 i;
		for(i=0;i<16;i++){Q[i]=0;}
		for(i=0;i<4;i++){R[i]=0;}

		for(i=0;i<4;i++){x_est[i]=0;}
		for(i=0;i<4;i++){x_pred[i]=0;}
		for(i=0;i<16;i++){P_pred[i]=0;}
		for(i=0;i<16;i++){P_est[i]=0;}
		P_est[0]=0x7FFF;
		P_est[5]=0x7FFF;
		P_est[10]=0x7FFF;
		P_est[15]=0x7FFF;
		for(i=0;i<2;i++){Y_mes[i]=0;}
		for(i=0;i<2;i++){ukalm[i]=0;}
		for(i=0;i<8;i++){Kalm[i]=0;}

		for(i=0;i<16;i++){PSI[i]=0;}

		set_dim(4);
		E._mu()=zeros(4);
		E._R()=zeros(4,4);
		init_ekf(0.000125);
	};
	//! Here dt = [yt;ut] of appropriate dimensions
	void bayes ( const vec &yt, const vec &ut );
	//!dummy!
	const epdf& posterior() const {return E;};
	
};

UIREGISTER(EKFfixed);

#endif

//! EKF for testing q44
class EKFtest: public EKF_UD{
	void bayes ( const vec &yt, const vec &cond ) {
		EKF_UD::bayes(yt,cond);
		vec D =  prior()._R()._D();
		
		if (D(3)>10) D(3) = 10;
		
		prior()._R().__D()=D;
	}
};
UIREGISTER(EKFtest);

/*!
\brief Extended Kalman Filter with UD matrices in fixed point16 arithmetic

An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean.
*/
class EKFfixedUD : public BM {
	public:
		LOG_LEVEL(EKFfixedUD,logU, logG, logD, logA, logP);
		
		void init_ekf(double Tv);
		void ekf(double ux, double uy, double isxd, double isyd);
				
		/* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/
		int16 Q[16]; /* matrix [4,4] */
		int16 R[4]; /* matrix [2,2] */
		
		int16 x_est[4]; /* estimate and prediction */
		
		int16 PSI[16]; /* matrix [4,4] */
		int16 PSIU[16]; /* matrix PIS*U, [4,4] */
		
		int16 Uf[16]; // upper triangular of covariance (inplace)
		int16 Df[4];  // diagonal covariance
		int16 Dfold[4]; // temp of D
		int16 G[16];  // temp for bierman
		
		int16 cA, cB, cC, cG, cH;  // cD, cE, cF, cI ... nepouzivane
		
		enorm<fsqmat> E;
		mat Ry;
		
	public:
		//! Default constructor
		EKFfixedUD ():BM(),E(),Ry(2,2){
			int16 i;
			for(i=0;i<16;i++){Q[i]=0;}
			for(i=0;i<4;i++){R[i]=0;}
			
			for(i=0;i<4;i++){x_est[i]=0;}
			for(i=0;i<16;i++){Uf[i]=0;}
			for(i=0;i<4;i++){Df[i]=0;}
			for(i=0;i<16;i++){G[i]=0;}
			for(i=0;i<4;i++){Dfold[i]=0;}
			
			for(i=0;i<16;i++){PSI[i]=0;}
			
			set_dim(4);
			dimy = 2;
			dimc = 2;
			E._mu()=zeros(4);
			E._R()=zeros(4,4);
			init_ekf(0.000125);
		};
		//! Here dt = [yt;ut] of appropriate dimensions
		void bayes ( const vec &yt, const vec &ut );
		//!dummy!
		const epdf& posterior() const {return E;};
		void log_register(logger &L, const string &prefix){
			BM::log_register ( L, prefix );
			
				L.add_vector ( log_level, logG, RV("G",16), prefix );
				L.add_vector ( log_level, logU, RV ("U", 16 ), prefix );
				L.add_vector ( log_level, logD, RV ("D", 4 ), prefix );
				L.add_vector ( log_level, logA, RV ("A", 16 ), prefix );
				L.add_vector ( log_level, logP, RV ("P", 16 ), prefix );
				
		};
		//void from_setting();
};

UIREGISTER(EKFfixedUD);

/*!
 * \brief Extended Kalman Filter with UD matrices in fixed point16 arithmetic
 * 
 * An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean.
 */
class EKFfixedUD2 : public BM {
public:
	LOG_LEVEL(EKFfixedUD2,logU, logG, logD, logA, logC, logP);
	
	void init_ekf2(double Tv);
	void ekf2(double ux, double uy, double isxd, double isyd);
	
	/* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/
	int16 Q[4]; /* matrix [4,4] */
	int16 R[4]; /* matrix [2,2] */
	
	int16 x_est[2]; /* estimate and prediction */
	int16 y_est[2]; /* estimate and prediction */
	int16 y_old[2]; /* estimate and prediction */
	
	int16 PSI[4]; /* matrix [4,4] */
	int16 PSIU[4]; /* matrix PIS*U, [4,4] */
	int16 C[4]; /* matrix [4,4] */
	
	int16 Uf[4]; // upper triangular of covariance (inplace)
	int16 Df[2];  // diagonal covariance
	int16 Dfold[2]; // temp of D
	int16 G[4];  // temp for bierman
	
	int16 cA, cB, cC, cG, cH;  // cD, cE, cF, cI ... nepouzivane
	
	enorm<fsqmat> E;
	mat Ry;
	
public:
	//! Default constructor
	EKFfixedUD2 ():BM(),E(),Ry(2,2){
		int16 i;
		for(i=0;i<4;i++){Q[i]=0;}
		for(i=0;i<4;i++){R[i]=0;}
		
		for(i=0;i<2;i++){x_est[i]=0;}
		for(i=0;i<2;i++){y_est[i]=0;}
		for(i=0;i<2;i++){y_old[i]=0;}
		for(i=0;i<4;i++){Uf[i]=0;}
		for(i=0;i<2;i++){Df[i]=0;}
		for(i=0;i<4;i++){G[i]=0;}
		for(i=0;i<2;i++){Dfold[i]=0;}
		
		for(i=0;i<4;i++){PSI[i]=0;}
		for(i=0;i<4;i++){C[i]=0;}
		
		set_dim(2);
		dimc = 2;
		dimy = 2;
		E._mu()=zeros(2);
		E._R()=zeros(2,2);
		init_ekf2(0.000125);
	};
	//! Here dt = [yt;ut] of appropriate dimensions
	void bayes ( const vec &yt, const vec &ut );
	//!dummy!
	const epdf& posterior() const {return E;};
	void log_register(logger &L, const string &prefix){
		BM::log_register ( L, prefix );
		
		L.add_vector ( log_level, logG, RV("G2",4), prefix );
		L.add_vector ( log_level, logU, RV ("U2", 4 ), prefix );
		L.add_vector ( log_level, logD, RV ("D2", 2 ), prefix );
		L.add_vector ( log_level, logA, RV ("A2", 4 ), prefix );
		L.add_vector ( log_level, logC, RV ("C2", 4 ), prefix );
		L.add_vector ( log_level, logP, RV ("P2", 4 ), prefix );
		
	};
	//void from_setting();
};

UIREGISTER(EKFfixedUD2);

/*!
 * \brief Extended Kalman Filter with UD matrices in fixed point16 arithmetic
 * 
 * An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean.
 */
class EKFfixedUD3 : public BM {
public:
	LOG_LEVEL(EKFfixedUD3,logU, logG, logD, logA, logC, logP);
	
	void init_ekf3(double Tv);
	void ekf3(double ux, double uy, double isxd, double isyd);
	
	/* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/
	int16 Q[9]; /* matrix [4,4] */
	int16 R[4]; /* matrix [2,2] */
	
	int16 x_est[3]; /* estimate and prediction */
	int16 y_est[2]; /* estimate and prediction */
	int16 y_old[2]; /* estimate and prediction */
	
	int16 PSI[9]; /* matrix [4,4] */
	int16 PSIU[9]; /* matrix PIS*U, [4,4] */
	int16 C[6]; /* matrix [4,4] */
	
	int16 Uf[9]; // upper triangular of covariance (inplace)
	int16 Df[3];  // diagonal covariance
	int16 Dfold[3]; // temp of D
	int16 G[9];  // temp for bierman
	
	int16 cA, cB, cC, cG, cF, cH;  // cD, cE, cF, cI ... nepouzivane
	
	enorm<fsqmat> E;
	mat Ry;
	
public:
	//! Default constructor
	EKFfixedUD3 ():BM(),E(),Ry(2,2){
		int16 i;
		for(i=0;i<9;i++){Q[i]=0;}
		for(i=0;i<4;i++){R[i]=0;}
		
		for(i=0;i<3;i++){x_est[i]=0;}
		for(i=0;i<2;i++){y_est[i]=0;}
		for(i=0;i<2;i++){y_old[i]=0;}
		for(i=0;i<9;i++){Uf[i]=0;}
		for(i=0;i<3;i++){Df[i]=0;}
		for(i=0;i<4;i++){G[i]=0;}
		for(i=0;i<3;i++){Dfold[i]=0;}
		
		for(i=0;i<9;i++){PSI[i]=0;}
		for(i=0;i<6;i++){C[i]=0;}
		
		set_dim(3);
		dimc = 2;
		dimy = 2;
		E._mu()=zeros(3);
		E._R()=zeros(3,3);
		init_ekf3(0.000125);
	};
	//! Here dt = [yt;ut] of appropriate dimensions
	void bayes ( const vec &yt, const vec &ut );
	//!dummy!
	const epdf& posterior() const {return E;};
	void log_register(logger &L, const string &prefix){
		BM::log_register ( L, prefix );		
	};
	//void from_setting();
};

UIREGISTER(EKFfixedUD3);

/*!
 * \brief Extended Kalman Filter with Chol matrices in fixed point16 arithmetic
 * 
 * An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean.
 */
class EKFfixedCh : public BM {
public:
	LOG_LEVEL(EKFfixedCh,logCh, logA, logP);
	
	void init_ekf(double Tv);
	void ekf(double ux, double uy, double isxd, double isyd);
	
	/* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/
	int16 Q[16]; /* matrix [4,4] */
	int16 R[4]; /* matrix [2,2] */
	
	int16 x_est[4]; /* estimate and prediction */
	
	int16 PSI[16]; /* matrix [4,4] */
	int16 PSICh[16]; /* matrix PIS*U, [4,4] */
	
	int16 Chf[16]; // upper triangular of covariance (inplace)
	
	int16 cA, cB, cC, cG, cH;  // cD, cE, cF, cI ... nepouzivane
	
	enorm<chmat> E;
	mat Ry;
	
public:
	//! Default constructor
	EKFfixedCh ():BM(),E(),Ry(2,2){
		int16 i;
		for(i=0;i<16;i++){Q[i]=0;}
		for(i=0;i<4;i++){R[i]=0;}
		
		for(i=0;i<4;i++){x_est[i]=0;}
		for(i=0;i<16;i++){Chf[i]=0;}
		
		for(i=0;i<16;i++){PSI[i]=0;}
		
		set_dim(4);
		dimc = 2;
		dimy =2;
		E._mu()=zeros(4);
		E._R()=zeros(4,4);
		init_ekf(0.000125);
	};
	//! Here dt = [yt;ut] of appropriate dimensions
	void bayes ( const vec &yt, const vec &ut );
	//!dummy!
	const epdf& posterior() const {return E;};
	void log_register(logger &L, const string &prefix){
		BM::log_register ( L, prefix );
		
		L.add_vector ( log_level, logCh, RV ("Ch", 16 ), prefix );
		L.add_vector ( log_level, logA, RV ("A", 16 ), prefix );
		L.add_vector ( log_level, logP, RV ("P", 16 ), prefix );
		
	};
	//void from_setting();
};

UIREGISTER(EKFfixedCh);

/*!
 * \brief Extended Kalman Filter with UD matrices in fixed point16 arithmetic
 * 
 * An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean.
 */
class EKFfixedCh2 : public BM {
public:
	LOG_LEVEL(EKFfixedCh2,logCh, logA, logC, logP);
	
	void init_ekf2(double Tv);
	void ekf2(double ux, double uy, double isxd, double isyd);
	
	/* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/
	int16 Q[4]; /* matrix [4,4] */
	int16 R[4]; /* matrix [2,2] */
	
	int16 x_est[2]; /* estimate and prediction */
	int16 y_est[2]; /* estimate and prediction */
	int16 y_old[2]; /* estimate and prediction */
	
	int16 PSI[4]; /* matrix [4,4] */
	int16 PSICh[4]; /* matrix PIS*U, [4,4] */
	int16 C[4]; /* matrix [4,4] */
	
	int16 Chf[4]; // upper triangular of covariance (inplace)
	
	int16 cA, cB, cC, cG, cH;  // cD, cE, cF, cI ... nepouzivane
	
	enorm<fsqmat> E;
	mat Ry;
	
public:
	//! Default constructor
	EKFfixedCh2 ():BM(),E(),Ry(2,2){
		int16 i;
		for(i=0;i<4;i++){Q[i]=0;}
		for(i=0;i<4;i++){R[i]=0;}
		
		for(i=0;i<2;i++){x_est[i]=0;}
		for(i=0;i<2;i++){y_est[i]=0;}
		for(i=0;i<2;i++){y_old[i]=0;}
		for(i=0;i<4;i++){Chf[i]=0;}
		
		for(i=0;i<4;i++){PSI[i]=0;}
		for(i=0;i<4;i++){C[i]=0;}
		
		set_dim(2);
		dimc = 2;
		dimy = 2;
		E._mu()=zeros(2);
		E._R()=zeros(2,2);
		init_ekf2(0.000125);
	};
	//! Here dt = [yt;ut] of appropriate dimensions
	void bayes ( const vec &yt, const vec &ut );
	//!dummy!
	const epdf& posterior() const {return E;};
	void log_register(logger &L, const string &prefix){
		BM::log_register ( L, prefix );
		
		L.add_vector ( log_level, logCh, RV ("Ch2", 4 ), prefix );
		L.add_vector ( log_level, logA, RV ("A2", 4 ), prefix );
		L.add_vector ( log_level, logC, RV ("C2", 4 ), prefix );
		L.add_vector ( log_level, logP, RV ("P2", 4 ), prefix );
		
	};
	//void from_setting();
};

UIREGISTER(EKFfixedCh2);


//! EKF for comparison of EKF_UD with its fixed-point16 implementation
class EKF_UDfix : public BM {
	protected:
		//! logger
		LOG_LEVEL(EKF_UDfix,logU, logG);
		//! Internal Model f(x,u)
		shared_ptr<diffbifn> pfxu;
		
		//! Observation Model h(x,u)
		shared_ptr<diffbifn> phxu;
		
		//! U part
		mat U;
		//! D part
		vec D;
				
		mat A;
		mat C;
		mat Q;
		vec R;
		
		enorm<ldmat> est;
		
		
	public:
		
		//! copy constructor duplicated 
		EKF_UDfix* _copy() const {
			return new EKF_UDfix(*this);
		}
		
		const enorm<ldmat>& posterior()const{return est;};
		
		enorm<ldmat>& prior() {
			return const_cast<enorm<ldmat>&>(posterior());
		}
		
		EKF_UDfix(){}
		
		
		EKF_UDfix(const EKF_UDfix &E0): pfxu(E0.pfxu),phxu(E0.phxu), U(E0.U), D(E0.D){}
		
		//! Set nonlinear functions for mean values and covariance matrices.
		void set_parameters ( const shared_ptr<diffbifn> &pfxu, const shared_ptr<diffbifn> &phxu, const mat Q0, const vec R0 );
		
		//! Here dt = [yt;ut] of appropriate dimensions
		void bayes ( const vec &yt, const vec &cond = empty_vec );
		
		void log_register ( bdm::logger& L, const string& prefix ){
			BM::log_register ( L, prefix );
			
			if ( log_level[logU] )
				L.add_vector ( log_level, logU, RV ( dimension()*dimension() ), prefix );
			if ( log_level[logG] )
				L.add_vector ( log_level, logG, RV ( dimension()*dimension() ), prefix );
			
		}
		/*! Create object from the following structure
		
		\code
		class = 'EKF_UD';
		OM = configuration of bdm::diffbifn;    % any offspring of diffbifn, bdm::diffbifn::from_setting
		IM = configuration of bdm::diffbifn;    % any offspring of diffbifn, bdm::diffbifn::from_setting
		dQ = [...];                             % vector containing diagonal of Q
		dR = [...];                             % vector containing diagonal of R
		--- optional fields ---
		mu0 = [...];                            % vector of statistics mu0
		dP0 = [...];                            % vector containing diagonal of P0
		-- or --
		P0 = [...];                             % full matrix P0
		--- inherited fields ---
		bdm::BM::from_setting
		\endcode
		If the optional fields are not given, they will be filled as follows:
		\code
		mu0 = [0,0,0,....];                     % empty statistics
		P0 = eye( dim );              
		\endcode
		*/
		void from_setting ( const Setting &set );
		
		void validate() {};
		// TODO dodelat void to_setting( Setting &set ) const;
		
};
UIREGISTER(EKF_UDfix);


class MPF_pmsm_red:public BM{
		double qom, qth, r;


public:
		MPF_pmsm_red(){
			dimy=2;
			dimc=2;
			qom=1e-1;
			qth=1e-6;
			r=1e-1;
		};
	void bayes ( const vec &val, const vec &cond ) {
		const double &isa = val(0);
		const double &isb = val(1);
		const double &usa = cond(0);
		const double &usb = cond(1);
		mpf_bayes(isa,isb,usa,usb);
	}
	
	class mp:public epdf{
	public:
		mp():epdf(){set_dim(3);}
		vec sample() const {return zeros(3);}
		double evallog(const vec &v) const {return 0.0;}
		vec mean() const {vec tmp(3); mpf_mean(tmp._data(), tmp._data()+1, tmp._data()+2); return tmp;}
		vec variance() const {return zeros(3);}
	};
	
	mp mypdf;
	const mp& posterior() const {return mypdf;}
	
	void from_setting(const Setting &set){
		BM::from_setting(set);
		UI::get(qom,set,"qom",UI::optional);
		UI::get(qth,set,"qth",UI::optional);
		UI::get(r,set,"r",UI::optional);
	}
	void validate(){
		mpf_init(qom,qth,r);
		
	}
};
UIREGISTER(MPF_pmsm_red);


#endif // KF_H


