#include "arx.h"
namespace bdm {

void ARX::bayes_weighted ( const vec &yt, const vec &cond, const double w ) {
	bdm_assert_debug ( yt.length() == dimy, "BM::bayes yt is of size "+num2str(yt.length())+" expected dimension is "+num2str(dimy) );
	bdm_assert_debug ( cond.length() == rgrlen , "BM::bayes cond is of size "+num2str(cond.length())+" expected dimension is "+num2str(rgrlen) );
	
	BMEF::bayes_weighted(yt,cond,w); //potential discount scheduling
	
	double lnc;
	//cache
	ldmat &V = est._V();
	double &nu = est._nu();

	dyad.set_subvector ( 0, yt );
	if (cond.length()>0)
		dyad.set_subvector ( dimy, cond );
	// possible "1" is there from the beginning

	if ( frg < 1.0 ) {
		est.pow ( frg ); // multiply V and nu


		//stabilize
		ldmat V0 = alter_est._V(); //$ copy
		double &nu0 = alter_est._nu();

		V0 *= ( 1 - frg );
		V += V0; //stabilization
		nu += ( 1 - frg ) * nu0;

		// recompute loglikelihood of new "prior"
		if ( evalll ) {
			last_lognc = est.lognc();
		}
	}
	V.opupdt ( dyad, 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 &yt, const vec &cond ) const {
	egiw pred ( est );
	ldmat &V = pred._V();
	double &nu = pred._nu();

	double lll;
	vec dyad_p = dyad;
	dyad_p.set_subvector ( 0, yt );
	dyad_p.set_subvector(dimy,cond);
	
	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 ( dyad_p, 1.0 );
	nu += 1.0;
	// log(sqrt(2*pi)) = 0.91893853320467
	return pred.lognc() - lll - 0.91893853320467;
}

void ARX::flatten ( const BMEF* B , double weight =1.0) {
	const ARX* A = dynamic_cast<const ARX*> ( B );
	// nu should be equal to B.nu
	est.pow ( A->posterior()._nu() / posterior()._nu() *weight);
	if ( evalll ) {
		last_lognc = est.lognc();
	}
}

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 ( dimension() == A0->dimension(), "Statistics of different dimensions" );
	set_statistics ( A0->dimensiony(), A0->posterior()._V(), A0->posterior()._nu() );
}

enorm<ldmat>* ARX::epredictor ( const vec &cond ) const {
	bdm_assert_debug ( cond.length() == rgrlen , "ARX::epredictor cond is of size "+num2str(cond.length())+" expected dimension is "+num2str(rgrlen) );
	
	mat mu ( dimy, posterior()._V().rows() - dimy );
	mat R ( dimy, dimy );

	vec ext_rgr;
	if (have_constant){
		ext_rgr = concat(cond,vec_1(1.0));
	} else {
		ext_rgr = cond;
	}
	
	enorm<ldmat>* tmp;
	tmp = new enorm<ldmat> ( );
	//TODO: too hackish
	if ( yrv._dsize() > 0 ) {
	}

	est.mean_mat ( mu, R ); //mu =
	//correction for student-t  -- TODO check if correct!!
	//R*=nu/(nu-2);
	if (mu.cols()>0) {// nonempty egiw
		mat p_mu = mu.T() * ext_rgr; 	//the result is one column
		tmp->set_parameters ( p_mu.get_col ( 0 ), ldmat ( R ) );
	} else {
		tmp->set_parameters ( zeros( R.rows() ), ldmat ( R ) );
	}
	if (dimy==yrv._dsize())
		tmp->set_rv(yrv);
	return tmp;
}

mlstudent* ARX::predictor_student ( ) const {
	const ldmat &V = posterior()._V();

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

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

	int end = V._L().rows() - 1;
	ldmat Lam ( V._L() ( dimy, end, dimy, end ), V._D() ( dimy, 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 ( dimy, dimc ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam );
		}
	} else {
		// no constant term
		tmp->set_parameters ( mu, zeros ( dimy ), 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 indices current indices
\return best likelihood in the structure below the given one
*/
double egiw_bestbelow ( egiw Eg, egiw Eg0, double Egll, ivec &indices ) { //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 tmpindices;
	ivec maxindices = indices;


	cout << "bb:(" << indices << ") 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
			tmpindices = indices;
			tmpindices.del ( i );
			//search for a better match in this substructure
			belll = egiw_bestbelow ( Eg, Eg0, tmpll, tmpindices );
			if ( belll > maxll ) { //better match found
				maxll = belll;
				maxindices = tmpindices;
			}
		}
	}
	indices = maxindices;
	return maxll;
}

ivec ARX::structure_est ( const 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 ( const egiw &est0 ) {
	//some stuff with beliefs etc.
	ivec belief = vec_1 ( 2 );        // default belief
	int nbest = 1;           // nbest: how many regressors are returned
	int nrep = 5;         // nrep: number of random repetions of structure estimation
	double lambda   = 0.9;
	int k = 2;
	
	Array<str_aux> o2;
	
	ivec ind = bdm::straux1(est._V(),est._nu(), est0._V(), est0._nu(), belief, nbest, nrep, lambda, k, o2);
	
	return ind;
}

void ARX::from_setting ( const Setting &set ) {
	BMEF::from_setting(set);
	
	UI::get (rgr, set, "rgr", UI::compulsory );
	
	dimy = yrv._dsize();
	bdm_assert(dimy>0,"ARX::yrv should not be empty");
	rgrlen = rgr._dsize();

	int constant;
	if ( !UI::get ( constant, set, "constant", UI::optional ) ) {
		have_constant = true;
	} else {
		have_constant = constant > 0;
	}
	dimc = rgrlen;
	rvc = rgr;

	//init
	shared_ptr<egiw> pri = UI::build<egiw> ( set, "prior", UI::optional );
	if (pri){
		set_prior(pri.get());
	} else {
		shared_ptr<egiw> post = UI::build<egiw> ( set, "posterior", UI::optional );
		set_prior(post.get());
	}
		
	
	shared_ptr<egiw> alt = UI::build<egiw> ( set, "alternative", UI::optional );
	if ( alt ) {
		bdm_assert ( alt->_dimx() == dimy, "alternative is not compatible" );
		bdm_assert ( alt->_V().rows() == dimy + rgrlen + int(have_constant==true), "alternative is not compatible" );
		alter_est.set_parameters ( alt->_dimx(), alt->_V(), alt->_nu() );
		alter_est.validate();
	} 
	// frg handled by BMEF

}

void ARX::set_prior(const epdf *pri){
	const egiw * eg=dynamic_cast<const egiw*>(pri);
	if ( eg ) {
		bdm_assert ( eg->_dimx() == dimy, "prior is not compatible" );
		bdm_assert ( eg->_V().rows() == dimy + rgrlen + int(have_constant==true), "prior is not compatible" );
		est.set_parameters ( eg->_dimx(), eg->_V(), eg->_nu() );
		est.validate();
	} else {
		est.set_parameters ( dimy, zeros ( dimy + rgrlen +int(have_constant==true)) );
		set_prior_default ( est );
	}
	//check alternative
	if (alter_est.dimension()!=dimension()){
		alter_est = est;
	}
}
}

