/*!
  \file
  \brief Bayesian Filtering for generalized autoregressive (ARX) model
  \author Vaclav Smidl.

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

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

#ifndef AR_H
#define AR_H

#include "../math/functions.h"
#include "../stat/exp_family.h"
#include "../base/user_info.h"
//#include "../estim/kalman.h"
#include "arx_straux.h"

namespace bdm {

/*!
* \brief Linear Autoregressive model with Gaussian noise

Regression of the following kind:
\f[
y_t = \theta_1 \psi_1 + \theta_2 + \psi_2 +\ldots + \theta_n \psi_n + r e_t
\f]
where unknown parameters \c rv are \f$[\theta r]\f$, regression vector \f$\psi=\psi(y_{1:t},u_{1:t})\f$ is a known function of past outputs and exogeneous variables \f$u_t\f$. Distrubances \f$e_t\f$ are supposed to be normally distributed:
\f[
e_t \sim \mathcal{N}(0,1).
\f]

See \ref tut_arx for mathematical treatment.

The easiest way how to use the class is:
\include arx_simple.cpp

	\todo sort out constant terms - bayes should accept vec without additional 1s
*/
class ARX: public BMEF {
protected:
	//! switch if constant is modelled or not
	bool have_constant;
	//! vector of dyadic update
	vec dyad;
	//! RV of regressor
	RV rgr;
	//! length of the regressor (without optional constant)
	int rgrlen;
	//! posterior density
	egiw est;
	//! Alternative estimate of parameters, used in stabilized forgetting, see [Kulhavy]
	egiw alter_est;
public:
	//! \name Constructors
	//!@{
	ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ),  have_constant ( true ), dyad(), rgrlen(),est(), alter_est() {};
	ARX ( const ARX &A0 ) : BMEF ( A0 ),  have_constant ( A0.have_constant ), dyad ( A0.dyad ),rgrlen(A0.rgrlen), est ( A0.est ), alter_est ( A0.alter_est ) { };

	ARX* _copy() const;

	void set_frg ( double frg0 ) {
		frg = frg0;
	}
	void set_constant ( bool const0 ) {
		have_constant = const0;
	}
	void set_statistics ( int dimy0, const ldmat V0, double nu0 = -1.0 ) {
		est.set_parameters ( dimy0, V0, nu0 );
		est.validate();
		last_lognc = est.lognc();
		dimy = dimy0;
	}
	//!@}

	//! Set sufficient statistics
	void set_statistics ( const BMEF* BM0 );

	//!\name Mathematical operations
	//!@{

	//! Weighted Bayes \f$ dt = [y_t psi_t] \f$.
	void bayes_weighted ( const vec &yt, const vec &cond = empty_vec, const double w = 1.0 );
	void bayes ( const vec &yt, const vec &cond = empty_vec ) {
		bayes_weighted ( yt, cond, 1.0 );
	};
	double logpred ( const vec &yt, const vec &cond ) const;
	void flatten ( const BMEF* B , double weight );
	//! Conditioned version of the predictor
	enorm<ldmat>* epredictor ( const vec &rgr ) const;
	//! conditional version of the predictor
	template<class sq_T>
	shared_ptr<mlnorm<sq_T> > ml_predictor() const;
	//! fast version of predicto
	template<class sq_T>
	void ml_predictor_update ( mlnorm<sq_T> &pred ) const;
	mlstudent* predictor_student() const;
	//! Brute force structure estimation.\return indices of accepted regressors.
	ivec structure_est ( const egiw &Eg0 );
	//! Smarter structure estimation by Ludvik Tesar.\return indices of accepted regressors.
	ivec structure_est_LT ( const egiw &Eg0 );
	//! reduce structure to the given ivec of matrix V
	void reduce_structure(ivec &inds_in_V){
		ldmat V = posterior()._V();
		if (max(inds_in_V)>=V.rows()) {bdm_error("Incompatible structure");}
		
		ldmat newV(V,inds_in_V);
		est.set_parameters(dimy,newV, posterior()._nu());
		
		if (have_constant){
			ivec rgr_elem= find(inds_in_V<(V.rows()-1)); // < -- find non-constant 
			rgr = rgr.subselect(rgr_elem);
			rgrlen = rgr_elem.length();
		} else{
			rgr = rgr.subselect(inds_in_V);
		}
		validate();
	}
	//!@}

	//!\name Access attributes
	//!@{
	//! return correctly typed posterior (covariant return)
	const egiw& posterior() const {
		return est;
	}
	//!@}

	/*! UI for ARX estimator

	\code
	class = 'ARX';
	yrv   = RV({names_of_dt} )                 // description of output variables
	rgr   = RV({names_of_regressors}, [-1,-2]} // description of regressor variables
	constant = 1;                              // 0/1 switch if the constant term is modelled or not

	--- optional ---
	prior = {class='egiw',...};                // Prior density, when given default is used instead
	alternative = {class='egiw',...};          // Alternative density in stabilized estimation, when not given prior is used

	frg = 1.0;                                 // forgetting, default frg=1.0

	rv  = RV({names_of_parameters}}            // description of parametetr names
											   // default: [""]
	\endcode
	*/
	void from_setting ( const Setting &set );

	void validate() {
		BMEF::validate();	
		est.validate();
		
		// When statistics is defined, it has priority
		if(posterior()._dimx()>0) {//statistics is assigned
			dimy = posterior()._dimx();
			rgrlen=posterior()._V().rows() - dimy - int ( have_constant == true );
			dimc = rgrlen;
		} else{  // statistics is not assigned - build it from dimy and rgrlen
			bdm_assert(dimy>0,"No way to validate egiw: empty statistics and empty dimy");
			est.set_parameters(dimy, zeros(dimy+rgrlen+int(have_constant==true)));
			set_prior_default(est);
		}
		if (alter_est.dimension()==0) alter_est=est;

		dyad = ones ( est._V().rows() );
	}
	//! function sets prior and alternative density
	void set_prior ( const epdf *prior );
	//! build default prior and alternative when all values are set
	void set_prior_default ( egiw &prior ) {
		//assume
		vec dV0 ( prior._V().rows() );
		dV0.set_subvector ( 0, prior._dimx() - 1, 1.0 );
		if (dV0.length()>prior._dimx())
			dV0.set_subvector ( prior._dimx(), dV0.length() - 1, 1e-5 );
		
		prior.set_parameters ( prior._dimx(), ldmat ( dV0 ) );
		prior.validate();
	}

	void to_setting ( Setting &set ) const
	{			
		BMEF::to_setting( set ); // takes care of rv, yrv, rvc
		UI::save(rgr, set, "rgr");
		int constant = have_constant ? 1 : 0;
		UI::save(constant, set, "constant");
		UI::save(&alter_est, set, "alternative");
		UI::save(&posterior(), set, "posterior");
		
	} 
	//! access function
	RV & _rgr() {return rgr;}
	bool _have_constant() {return have_constant;}
	int _rgrlen() {return rgrlen;}
};
UIREGISTER ( ARX );
SHAREDPTR ( ARX );

/*! ARX model conditined by knowledge of the forgetting factor
\f[ f(\theta| d_1 \ldots d_t , \phi_t) \f]

The symbol \f$ \phi \f$ is assumed to be the last of the conditioning variables.
*/
class ARXfrg : public ARX {
public:
	ARXfrg() : ARX() {};
	//! copy constructor
	ARXfrg ( const ARXfrg &A0 ) : ARX ( A0 ) {};
	virtual ARXfrg* _copy() const {
		ARXfrg *A = new ARXfrg ( *this );
		return A;
	}

	void bayes ( const vec &val, const vec &cond ) {
		bdm_assert_debug(cond.size()>rgrlen, "ARXfrg: Insufficient conditioning, frg not given.");
		frg = cond ( rgrlen); // the first part after rgrlen
		ARX::bayes ( val, cond.left(rgrlen) );
	}
	void validate() {
		ARX::validate();
		rvc.add ( RV ( "{phi }", vec_1 ( 1 ) ) );
		dimc += 1;
	}
};
UIREGISTER ( ARXfrg );


/*!
* \brief ARX with partial forgetting

Implements partial forgetting when <tt>bayes(const vec &yt, const vec &cond=empty_vec)</tt> is called, where \c cond is a vector <em>(regressor', forg.factor')</em>.

Note, that the weights have the same order as hypotheses, following this scheme:
\li 0 means that the parameter doesn't change,
\li 1 means that the parameter varies.

The combinations of parameters are described binary:
\f[
  0, 0, 0\ldots \\
  1, 0, 0\ldots \\
  0, 1, 0\ldots \\
  1, 1, 0\ldots \\
  \ldots
\f]
where each n-th column has altering 1's and 0's in n-tuples, n = 0,...,number of params. Hence, the first forg. factor relates to the situation when no parameter changes, the second one when the first parameter changes etc.

See ARX class for more information about ARX.
*/
class ARXpartialforg : public ARX {
public:
	ARXpartialforg() : ARX(1.0) {};
	//! copy constructor
	ARXpartialforg ( const ARXpartialforg &A0 ) : ARX ( A0 ) {};
	virtual ARXpartialforg* _copy() const {
		ARXpartialforg *A = new ARXpartialforg ( *this );
		return A;
	}

    void bayes ( const vec &val, const vec &cond ); 
	void validate() {
		ARX::validate();
		int philen = 1 << (est._V().cols() - 1);
		rvc.add ( RV ( "{phi }", vec_1(philen) ) ); // pocet 2^parametru
		dimc += philen;
	}
};
UIREGISTER ( ARXpartialforg );


////////////////////
template<class sq_T>
shared_ptr< mlnorm<sq_T> > ARX::ml_predictor ( ) const {
	shared_ptr< mlnorm<sq_T> > tmp = new mlnorm<sq_T> ( );
	tmp->set_rv ( yrv );
	tmp->set_rvc ( _rvc() );

	ml_predictor_update ( *tmp );
	tmp->validate();
	return tmp;
}

template<class sq_T>
void ARX::ml_predictor_update ( mlnorm<sq_T> &pred ) const {
	mat mu ( dimy, posterior()._V().rows() - dimy );
	mat R ( dimy, dimy );

	est.mean_mat ( mu, R ); //mu =
	mu = mu.T();
	//correction for student-t  -- TODO check if correct!!

	if ( have_constant ) { // constant term
		//Assume the constant term is the last one:
		pred.set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), sq_T ( R ) );
	} else {
		pred.set_parameters ( mu, zeros ( dimy ), sq_T ( R ) );
	}
}

};
#endif // AR_H