#include "arx.h"
namespace bdm {

void ARX::bayes ( const vec &dt, const double w ) {
	double lnc;
		
	if ( frg < 1.0 ) {
		est.pow ( frg );
		
		//stabilize
		ldmat V0(eye(V.rows()));
		V0*=(1-frg)*1e-3;
		V += V0; //stabilization
		nu +=(1-frg)*(0.1 + V.rows() + 1* dimx + 2);
		
		// recompute loglikelihood of "prior"
		if ( evalll ) {
			last_lognc = est.lognc();
		}
	}
	if (have_constant) {
		_dt.set_subvector(0,dt);
		V.opupdt ( _dt, w );
	} else {
		V.opupdt ( dt, w );
	}
	nu += w;

	// log(sqrt(2*pi)) = 0.91893853320467
	if ( evalll ) {
		lnc = est.lognc();
		ll = lnc - last_lognc - 0.91893853320467;
		last_lognc = lnc;
	}
}

double ARX::logpred ( const vec &dt ) const {
	egiw pred ( est );
	ldmat &V = pred._V();
	double &nu = pred._nu();

	double lll;

	if ( frg < 1.0 ) {
		pred.pow ( frg );
		lll = pred.lognc();
	} else//should be save: last_lognc is changed only by bayes;
		if ( evalll ) {
			lll = last_lognc;
		} else {
			lll = pred.lognc();
		}

	V.opupdt ( dt, 1.0 );
	nu += 1.0;
	// log(sqrt(2*pi)) = 0.91893853320467
	return pred.lognc() - lll - 0.91893853320467;
}

ARX* ARX::_copy_ ( ) const {
	ARX* Tmp = new ARX ( *this );
	return Tmp;
}

void ARX::set_statistics ( const BMEF* B0 ) {
	const ARX* A0 = dynamic_cast<const ARX*> ( B0 );

	bdm_assert_debug ( V.rows() == A0->V.rows(), "ARX::set_statistics Statistics  differ" );
	set_statistics ( A0->dimx, A0->V, A0->nu );
}

enorm<ldmat>* ARX::epredictor ( const vec &rgr ) const {
	int dim = dimx;//est.dimension();
	mat mu ( dim, V.rows() - dim );
	mat R ( dim, dim );

	enorm<ldmat>* tmp;
	tmp = new enorm<ldmat> ( );
	//TODO: too hackish
	if ( drv._dsize() > 0 ) {
	}

	est.mean_mat ( mu, R ); //mu =
	//correction for student-t  -- TODO check if correct!!
	//R*=nu/(nu-2);
	mat p_mu = mu.T() * rgr; 	//the result is one column
	tmp->set_parameters ( p_mu.get_col ( 0 ), ldmat ( R ) );
	return tmp;
}

mlnorm<ldmat>* ARX::predictor ( ) const {
	int dim = est.dimension();
	
	mat mu ( dim, V.rows() - dim );
	mat R ( dim, dim );
	mlnorm<ldmat>* tmp;
	tmp = new mlnorm<ldmat> ( );

	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:
		tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ) );
	} else {
		tmp->set_parameters ( mu, zeros ( dim ), ldmat ( R ) );
	}
	return tmp;
}

mlstudent* ARX::predictor_student ( ) const {
	int dim = est.dimension();

	mat mu ( dim, V.rows() - dim );
	mat R ( dim, dim );
	mlstudent* tmp;
	tmp = new mlstudent ( );

	est.mean_mat ( mu, R ); //
	mu = mu.T();

	int xdim = dimx;
	int end = V._L().rows() - 1;
	ldmat Lam ( V._L() ( xdim, end, xdim, end ), V._D() ( xdim, end ) );  //exp val of R


	if ( have_constant) { // no constant term
		//Assume the constant term is the last one:
		if ( mu.cols() > 1 ) {
			tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam );
		} else {
			tmp->set_parameters ( mat ( dim, 0 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam );
		}
	} else {
		// no constant term
		tmp->set_parameters ( mu, zeros ( xdim ), ldmat ( R ), Lam );
	}
	return tmp;
}



/*! \brief Return the best structure
@param Eg a copy of GiW density that is being examined
@param Eg0 a copy of prior GiW density before estimation
@param Egll likelihood of the current Eg
@param indeces current indeces
\return best likelihood in the structure below the given one
*/
double egiw_bestbelow ( egiw Eg, egiw Eg0, double Egll, ivec &indeces ) { //parameter Eg is a copy!
	ldmat Vo = Eg._V(); //copy
	ldmat Vo0 = Eg._V(); //copy
	ldmat& Vp = Eg._V(); // pointer into Eg
	ldmat& Vp0 = Eg._V(); // pointer into Eg
	int end = Vp.rows() - 1;
	int i;
	mat Li;
	mat Li0;
	double maxll = Egll;
	double tmpll = Egll;
	double belll = Egll;

	ivec tmpindeces;
	ivec maxindeces = indeces;


	cout << "bb:(" << indeces << ") ll=" << Egll << endl;

	//try to remove only one rv
	for ( i = 0; i < end; i++ ) {
		//copy original
		Li = Vo._L();
		Li0 = Vo0._L();
		//remove stuff
		Li.del_col ( i + 1 );
		Li0.del_col ( i + 1 );
		Vp.ldform ( Li, Vo._D() );
		Vp0.ldform ( Li0, Vo0._D() );
		tmpll = Eg.lognc() - Eg0.lognc(); // likelihood is difference of norm. coefs.

		cout << "i=(" << i << ") ll=" << tmpll << endl;

		//
		if ( tmpll > Egll ) { //increase of the likelihood
			tmpindeces = indeces;
			tmpindeces.del ( i );
			//search for a better match in this substructure
			belll = egiw_bestbelow ( Eg, Eg0, tmpll, tmpindeces );
			if ( belll > maxll ) { //better match found
				maxll = belll;
				maxindeces = tmpindeces;
			}
		}
	}
	indeces = maxindeces;
	return maxll;
}

ivec ARX::structure_est ( egiw est0 ) {
	ivec ind = linspace ( 1, est.dimension() - 1 );
	egiw_bestbelow ( est, est0, est.lognc() - est0.lognc(), ind );
	return ind;
}



ivec ARX::structure_est_LT ( egiw est0 ) {
	//some stuff with beliefs etc.
//	ivec ind = bdm::straux1(V,nu, est0._V(), est0._nu());
	return ivec();//ind;
}

void ARX::from_setting ( const Setting &set ) {
	shared_ptr<RV> yrv = UI::build<RV> ( set, "rv", UI::compulsory );
	shared_ptr<RV> rrv = UI::build<RV> ( set, "rgr", UI::compulsory );
	int ylen = yrv->_dsize();
	// rgrlen - including constant!!!
	int rgrlen = rrv->_dsize();
	
	set_rv ( *yrv, *rrv );
	
	string opt;
	if ( UI::get(opt, set,  "options", UI::optional) ) {
		BM::set_options(opt);
	}
	int constant;
	if (!UI::get(constant, set, "constant", UI::optional)){
		have_constant=true;
	} else {
		have_constant=constant>0;
	}
	if (have_constant) {rgrlen++;_dt=ones(rgrlen+ylen);}

	//init
	mat V0;
	vec dV0;
	if (!UI::get(V0, set, "V0",UI::optional)){
		if ( !UI::get ( dV0, set, "dV0" ) )
			dV0 = concat ( 1e-3 * ones ( ylen ), 1e-5 * ones ( rgrlen ) );
		V0 = diag ( dV0 );
	}
	double nu0;
	if ( !UI::get ( nu0, set, "nu0" ) )
		nu0 = rgrlen + ylen + 2;

	double frg;
	if ( !UI::get ( frg, set, "frg" ) )
		frg = 1.0;

	set_parameters ( frg );
	set_statistics ( ylen, V0, nu0 );
	
	//name results (for logging)
	shared_ptr<RV> rv_par=UI::build<RV>(set, "rv_param",UI::optional );
	if (!rv_par){
		est.set_rv ( RV ( "{theta r }", vec_2 ( ylen*rgrlen, ylen*ylen ) ) );
	} else {
		est.set_rv ( *rv_par );
	}
	validate();
}

}