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

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

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

#ifndef KF_H
#define KF_H


#include "../math/functions.h"
#include "../stat/exp_family.h"
#include "../math/chmat.h"
#include "../base/user_info.h"

namespace bdm {

/*!
* \brief Basic Kalman filter with full matrices (education purpose only)! Will be deleted soon!
*/

class KalmanFull {
protected:
	int dimx, dimy, dimu;
	mat A, B, C, D, R, Q;

	//cache
	mat _Pp, _Ry, _iRy, _K;
public:
	//posterior
	//! Mean value of the posterior density
	vec mu;
	//! Variance of the posterior density
	mat P;

	bool evalll;
	double ll;
public:
	//! Full constructor
	KalmanFull ( mat A, mat B, mat C, mat D, mat R, mat Q, mat P0, vec mu0 );
	//! Here dt = [yt;ut] of appropriate dimensions
	void bayes ( const vec &dt );
	//! print elements of KF
	friend std::ostream &operator<< ( std::ostream &os, const KalmanFull &kf );
	//! For EKFfull;
	KalmanFull() {};
};


/*!
* \brief Kalman filter with covariance matrices in square root form.

Parameter evolution model:\f[ x_t = A x_{t-1} + B u_t + Q^{1/2} e_t \f]
Observation model: \f[ y_t = C x_{t-1} + C u_t + Q^{1/2} w_t. \f]
Where $e_t$ and $w_t$ are independent vectors Normal(0,1)-distributed disturbances.
*/
template<class sq_T>

class Kalman : public BM {
protected:
	//! Indetifier of output rv
	RV rvy;
	//! Indetifier of exogeneous rv
	RV rvu;
	//! cache of rv.count()
	int dimx;
	//! cache of rvy.count()
	int dimy;
	//! cache of rvu.count()
	int dimu;
	//! Matrix A
	mat A;
	//! Matrix B
	mat B;
	//! Matrix C
	mat C;
	//! Matrix D
	mat D;
	//! Matrix Q in square-root form
	sq_T Q;
	//! Matrix R in square-root form
	sq_T R;

	//!posterior density on $x_t$
	enorm<sq_T> est;
	//!preditive density on $y_t$
	enorm<sq_T> fy;

	//! placeholder for Kalman gain
	mat _K;
	//! cache of fy.mu
	vec& _yp;
	//! cache of fy.R
	sq_T& _Ry;
	//!cache of est.mu
	vec& _mu;
	//!cache of est.R
	sq_T& _P;

public:
	//! Default constructor
	Kalman ( );
	//! Copy constructor
	Kalman ( const Kalman<sq_T> &K0 );
	//! Set parameters with check of relevance
	void set_parameters ( const mat &A0, const mat &B0, const mat &C0, const mat &D0, const sq_T &Q0, const sq_T &R0 );
	//! Set estimate values, used e.g. in initialization.
	void set_est ( const vec &mu0, const sq_T &P0 ) {
		sq_T pom ( dimy );
		est.set_parameters ( mu0, P0 );
		P0.mult_sym ( C, pom );
		fy.set_parameters ( C*mu0, pom );
	};

	//! Here dt = [yt;ut] of appropriate dimensions
	void bayes ( const vec &dt );
	//!access function
	const epdf& posterior() const {
		return est;
	}
	//!access function
	mat& __K() {
		return _K;
	}
	//!access function
	vec _dP() {
		return _P->getD();
	}
};

/*! \brief Kalman filter in square root form

Trivial example:
\include kalman_simple.cpp

*/
class KalmanCh : public Kalman<chmat> {
protected:
//! pre array (triangular matrix)
	mat preA;
//! post array (triangular matrix)
	mat postA;

public:
	//! copy constructor
	BM* _copy_() const {
		KalmanCh* K = new KalmanCh;
		K->set_parameters ( A, B, C, D, Q, R );
		K->set_statistics ( _mu, _P );
		return K;
	}
	//! Set parameters with check of relevance
	void set_parameters ( const mat &A0, const mat &B0, const mat &C0, const mat &D0, const chmat &Q0, const chmat &R0 );
	void set_statistics ( const vec &mu0, const chmat &P0 ) {
		est.set_parameters ( mu0, P0 );
	};


	/*!\brief  Here dt = [yt;ut] of appropriate dimensions

	The following equality hold::\f[
	\left[\begin{array}{cc}
	R^{0.5}\\
	P_{t|t-1}^{0.5}C' & P_{t|t-1}^{0.5}CA'\\
	& Q^{0.5}\end{array}\right]<\mathrm{orth.oper.}>=\left[\begin{array}{cc}
	R_{y}^{0.5} & KA'\\
	& P_{t+1|t}^{0.5}\\
	\\\end{array}\right]\f]

	Thus this object evaluates only predictors! Not filtering densities.
	*/
	void bayes ( const vec &dt );
};

/*!
\brief Extended Kalman Filter in full matrices

An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean.
*/
class EKFfull : public KalmanFull, public BM {
protected:
	//! Internal Model f(x,u)
	shared_ptr<diffbifn> pfxu;

	//! Observation Model h(x,u)
	shared_ptr<diffbifn> phxu;

	enorm<fsqmat> E;
public:
	//! Default constructor
	EKFfull ( );

	//! 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 mat R0 );

	//! Here dt = [yt;ut] of appropriate dimensions
	void bayes ( const vec &dt );
	//! set estimates
	void set_statistics ( vec mu0, mat P0 ) {
		mu = mu0;
		P = P0;
	};
	//!dummy!
	const epdf& posterior() const {
		return E;
	};
	const mat _R() {
		return P;
	}
};

/*!
\brief Extended Kalman Filter

An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean.
*/
template<class sq_T>
class EKF : public Kalman<fsqmat> {
	//! Internal Model f(x,u)
	diffbifn* pfxu;
	//! Observation Model h(x,u)
	diffbifn* phxu;
public:
	//! Default constructor
	EKF ( RV rvx, RV rvy, RV rvu );
	//! copy constructor
	EKF<sq_T>* _copy_() const {
		return new EKF<sq_T> ( this );
	}
	//! Set nonlinear functions for mean values and covariance matrices.
	void set_parameters ( diffbifn* pfxu, diffbifn* phxu, const sq_T Q0, const sq_T R0 );
	//! Here dt = [yt;ut] of appropriate dimensions
	void bayes ( const vec &dt );
};

/*!
\brief Extended Kalman Filter in Square root

An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean.
*/

class EKFCh : public KalmanCh {
protected:
	//! Internal Model f(x,u)
	shared_ptr<diffbifn> pfxu;

	//! Observation Model h(x,u)
	shared_ptr<diffbifn> phxu;
public:
	//! copy constructor duplicated - calls different set_parameters
	BM* _copy_() const {
		EKFCh* E = new EKFCh;
		E->set_parameters ( pfxu, phxu, Q, R );
		E->set_statistics ( _mu, _P );
		return E;
	}
	//! Set nonlinear functions for mean values and covariance matrices.
	void set_parameters ( const shared_ptr<diffbifn> &pfxu, const shared_ptr<diffbifn> &phxu, const chmat Q0, const chmat R0 );

	//! Here dt = [yt;ut] of appropriate dimensions
	void bayes ( const vec &dt );

	void from_setting ( const Setting &set );

	// TODO dodelat void to_setting( Setting &set ) const;

	const enorm<chmat>& posterior() {return est;}
};

UIREGISTER ( EKFCh );
SHAREDPTR ( EKFCh );


/*!
\brief Kalman Filter with conditional diagonal matrices R and Q.
*/

class KFcondQR : public Kalman<ldmat> {
//protected:
public:
	void condition ( const vec &QR ) {
		it_assert_debug ( QR.length() == ( dimx + dimy ), "KFcondRQ: conditioning by incompatible vector" );

		Q.setD ( QR ( 0, dimx - 1 ) );
		R.setD ( QR ( dimx, -1 ) );
	};
};

/*!
\brief Kalman Filter with conditional diagonal matrices R and Q.
*/

class KFcondR : public Kalman<ldmat> {
//protected:
public:
	//!Default constructor
	KFcondR ( ) : Kalman<ldmat> ( ) {};

	void condition ( const vec &R0 ) {
		it_assert_debug ( R0.length() == ( dimy ), "KFcondR: conditioning by incompatible vector" );

		R.setD ( R0 );
	};

};

//////// INstance

template<class sq_T>
Kalman<sq_T>::Kalman ( const Kalman<sq_T> &K0 ) : BM ( ), rvy ( K0.rvy ), rvu ( K0.rvu ),
		dimx ( K0.dimx ), dimy ( K0.dimy ), dimu ( K0.dimu ),
		A ( K0.A ), B ( K0.B ), C ( K0.C ), D ( K0.D ),
		Q ( K0.Q ), R ( K0.R ),
		est ( K0.est ), fy ( K0.fy ), _yp ( fy._mu() ), _Ry ( fy._R() ), _mu ( est._mu() ), _P ( est._R() ) {

// copy values in pointers
// 	_mu = K0._mu;
// 	_P = K0._P;
// 	_yp = K0._yp;
// 	_Ry = K0._Ry;

}

template<class sq_T>
Kalman<sq_T>::Kalman ( ) : BM (), est ( ), fy (),  _yp ( fy._mu() ), _Ry ( fy._R() ), _mu ( est._mu() ), _P ( est._R() ) {
};

template<class sq_T>
void Kalman<sq_T>::set_parameters ( const mat &A0, const  mat &B0, const mat &C0, const mat &D0, const sq_T &Q0, const sq_T &R0 ) {
	dimx = A0.rows();
	dimy = C0.rows();
	dimu = B0.cols();

	it_assert_debug ( A0.cols() == dimx, "Kalman: A is not square" );
	it_assert_debug ( B0.rows() == dimx, "Kalman: B is not compatible" );
	it_assert_debug ( C0.cols() == dimx, "Kalman: C is not square" );
	it_assert_debug ( ( D0.rows() == dimy ) || ( D0.cols() == dimu ),	"Kalman: D is not compatible" );
	it_assert_debug ( ( R0.cols() == dimy ) || ( R0.rows() == dimy ), "Kalman: R is not compatible" );
	it_assert_debug ( ( Q0.cols() == dimx ) || ( Q0.rows() == dimx ), "Kalman: Q is not compatible" );

	A = A0;
	B = B0;
	C = C0;
	D = D0;
	R = R0;
	Q = Q0;
}

template<class sq_T>
void Kalman<sq_T>::bayes ( const vec &dt ) {
	it_assert_debug ( dt.length() == ( dimy + dimu ), "KalmanFull::bayes wrong size of dt" );

	sq_T iRy ( dimy );
	vec u = dt.get ( dimy, dimy + dimu - 1 );
	vec y = dt.get ( 0, dimy - 1 );
	//Time update
	_mu = A * _mu + B * u;
	//P  = A*P*A.transpose() + Q; in sq_T
	_P.mult_sym ( A );
	_P  += Q;

	//Data update
	//_Ry = C*P*C.transpose() + R; in sq_T
	_P.mult_sym ( C, _Ry );
	_Ry  += R;

	mat Pfull = _P.to_mat();

	_Ry.inv ( iRy ); // result is in _iRy;
	_K = Pfull * C.transpose() * ( iRy.to_mat() );

	sq_T pom ( ( int ) Pfull.rows() );
	iRy.mult_sym_t ( C*Pfull, pom );
	( _P ) -= pom; // P = P -PC'iRy*CP;
	( _yp ) = C * _mu  + D * u; //y prediction
	( _mu ) += _K * ( y - _yp );


	if ( evalll == true ) { //likelihood of observation y
		ll = fy.evallog ( y );
	}

//cout << "y: " << y-(*_yp) <<" R: " << _Ry->to_mat() << " iR: " << _iRy->to_mat() << " ll: " << ll <<endl;

};

/*! \brief (Switching) Multiple Model
The model runs several models in parallel and evaluates thier weights (fittness).

The statistics of the resulting density are merged using (geometric?) combination.

The next step is performed with the new statistics for all models.
*/
class MultiModel: public BM {
protected:
	//! List of models between which we switch
	Array<EKFCh*> Models;
	//! vector of model weights
	vec w;
	//! cache of model lls
	vec _lls;
	//! type of switching policy [1=maximum,2=...]
	int policy;
	//! internal statistics
	enorm<chmat> est;
public:
	void set_parameters ( Array<EKFCh*> A, int pol0 = 1 ) {
		Models = A;//TODO: test if evalll is set
		w.set_length ( A.length() );
		_lls.set_length ( A.length() );
		policy = pol0;

		est.set_rv ( RV ( "MM", A ( 0 )->posterior().dimension(), 0 ) );
		est.set_parameters ( A ( 0 )->posterior().mean(), A ( 0 )->posterior()._R() );
	}
	void bayes ( const vec &dt ) {
		int n = Models.length();
		int i;
		for ( i = 0; i < n; i++ ) {
			Models ( i )->bayes ( dt );
			_lls ( i ) = Models ( i )->_ll();
		}
		double mlls = max ( _lls );
		w = exp ( _lls - mlls );
		w /= sum ( w ); //normalization
		//set statistics
		switch ( policy ) {
		case 1: {
			int mi = max_index ( w );
			const enorm<chmat> &st = Models ( mi )->posterior() ;
			est.set_parameters ( st.mean(), st._R() );
		}
		break;
		default:
			it_error ( "unknown policy" );
		}
		// copy result to all models
		for ( i = 0; i < n; i++ ) {
			Models ( i )->set_statistics ( est.mean(), est._R() );
		}
	}
	//! posterior density
	const enorm<chmat>& posterior() const {
		return est;
	}

	void from_setting ( const Setting &set );

	// TODO dodelat void to_setting( Setting &set ) const;

};

UIREGISTER ( MultiModel );
SHAREDPTR ( MultiModel );



//TODO why not const pointer??

template<class sq_T>
EKF<sq_T>::EKF ( RV rvx0, RV rvy0, RV rvu0 ) : Kalman<sq_T> ( rvx0, rvy0, rvu0 ) {}

template<class sq_T>
void EKF<sq_T>::set_parameters ( diffbifn* pfxu0,  diffbifn* phxu0, const sq_T Q0, const sq_T R0 ) {
	pfxu = pfxu0;
	phxu = phxu0;

	//initialize matrices A C, later, these will be only updated!
	pfxu->dfdx_cond ( _mu, zeros ( dimu ), A, true );
//	pfxu->dfdu_cond ( *_mu,zeros ( dimu ),B,true );
	B.clear();
	phxu->dfdx_cond ( _mu, zeros ( dimu ), C, true );
//	phxu->dfdu_cond ( *_mu,zeros ( dimu ),D,true );
	D.clear();

	R = R0;
	Q = Q0;
}

template<class sq_T>
void EKF<sq_T>::bayes ( const vec &dt ) {
	it_assert_debug ( dt.length() == ( dimy + dimu ), "KalmanFull::bayes wrong size of dt" );

	sq_T iRy ( dimy, dimy );
	vec u = dt.get ( dimy, dimy + dimu - 1 );
	vec y = dt.get ( 0, dimy - 1 );
	//Time update
	_mu = pfxu->eval ( _mu, u );
	pfxu->dfdx_cond ( _mu, u, A, false ); //update A by a derivative of fx

	//P  = A*P*A.transpose() + Q; in sq_T
	_P.mult_sym ( A );
	_P += Q;

	//Data update
	phxu->dfdx_cond ( _mu, u, C, false ); //update C by a derivative hx
	//_Ry = C*P*C.transpose() + R; in sq_T
	_P.mult_sym ( C, _Ry );
	( _Ry ) += R;

	mat Pfull = _P.to_mat();

	_Ry.inv ( iRy ); // result is in _iRy;
	_K = Pfull * C.transpose() * ( iRy.to_mat() );

	sq_T pom ( ( int ) Pfull.rows() );
	iRy.mult_sym_t ( C*Pfull, pom );
	( _P ) -= pom; // P = P -PC'iRy*CP;
	_yp = phxu->eval ( _mu, u ); //y prediction
	( _mu ) += _K * ( y - _yp );

	if ( evalll == true ) {
		ll += fy.evallog ( y );
	}
};


}
#endif // KF_H


