#include <itpp/itcomm.h>
#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!=1 ) mexErrMsgTxt ( "Wrong number of input variables!" );

	// Convert input variables to IT++ format
	mat Data = mxArray2mat ( input[0] );

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

	mat V0=1e-8*eye ( npsi ); V0 ( 0,0 ) *=100.0;
	double nu0=npsi;
	ARX Ar ( thr,V0,nu0 );
	// estimate
	for (int i=0;i<ndat;i++ ) {Ar.bayes ( Data.get_col ( i ) );	}

	vec th = Ar._epdf().mean();

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

	// Create output vectors
	output[0] = mxCreateDoubleMatrix ( 1,th.length(), mxREAL );
	// Convert the IT++ format to Matlab format for output
	vec2mxArray ( th, output[0] );

	if ( n_output>1 ) {
		double ll=Ar._tll();
		output[1] = mxCreateDoubleMatrix ( 1,1, mxREAL );
		double2mxArray ( ll, output[1] );
	}
	
	if ( n_output>2 ) {
		ivec str = Ar.structure_est ( egiw ( thr,V0,nu0 ) );
		output[2] = mxCreateDoubleMatrix ( 1,str.length(), mxREAL );
		ivec2mxArray ( str, output[2] );
	}
}