#include <itpp/itmex.h>

#include "../../bdm/estim/arx.h"

using namespace itpp;

void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) {
	// Check the number of inputs and output arguments
	if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output variables!" );
	if ( n_input!=3 ) mexErrMsgTxt ( "Wrong number of input variables!" );

	// Convert input variables to IT++ format
	mat Data = mxArray2mat ( input[0] );
	if (Data.rows()!=2)
		{ if (Data.cols()==2) Data=Data.T(); else mexErrMsgTxt("Data are not 2D!");}
	int ilow = mxArray2int ( input[1] )-1;
	int ihi = mxArray2int ( input[2] )-1; //coreection for different indeces

	// ------------------ Start of routine ---------------------------
	// ------------------ Start of routine ---------------------------
	int ndat=Data.cols();
	int npsi=Data.rows()+1; //add a constant
	RV thr ( "1","{ theta_r }",vec_1 ( npsi ),vec_1 ( 0 ) );

	mat V0=1e-8*eye ( npsi ); V0 ( 0,0 ) =10.0;
	vec rgr(npsi);
	double nu0=2;
	// fitting a linear part => third coef is "1"
	rgr(2) = 1.0;

	// RESULTS
	mat Tlls=zeros(ihi,ihi);
	// AR model
	ARX Ar ( thr,V0,nu0 );

	//Initialize
	for (int i=ilow; i<ihi; i++) {
		rgr.set_subvector(0,Data.get_col ( i ));
		Ar.bayes(rgr);
	}
	Ar.get_parameters(V0,nu0);
	V0 /=V0(1,1);
	nu0=3;
	Ar.set_parameters(V0,nu0);
		
	// FLATTEN
	for (int i=ilow; i<ihi; i++) {
		Ar.set_parameters(V0,nu0);
		for (int t=i+1; t<ihi; t++) {
			rgr.set_subvector(0,Data.get_col ( t ));
			Ar.bayes(rgr);
			Tlls(i,t)=Ar._tll();
		}
	}

	max_index(Tlls,ilow,ihi);

	// ------------------ End of routine -----------------------------

	output[0] = mxCreateDoubleMatrix ( 2,1, mxREAL );
	ivec2mxArray (vec_2(ilow+1,ihi+1) , output[0] );

	if ( n_output>1 ) {
	// Create output vectors
	output[1] = mxCreateDoubleMatrix ( Tlls.rows(),Tlls.cols(), mxREAL );
	// Convert the IT++ format to Matlab format for output
	mat2mxArray ( Tlls, output[1] );
	}
	
}