root/matlab/mex/linefit2.cpp @ 248

#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 ) =1e-4;
        vec rgr ( npsi );
        double nu0=0.1;
        // fitting a linear part => third coef is "1"
        rgr ( 2 ) = 1.0;

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

        //Initialize
        for ( int i=ilow+2; i<ihi-2; i++ ) {
                rgr.set_subvector ( 0,Data.get_col ( i ) );
                Ar.bayes ( rgr );
        }
        // FLATTEN
/*      Ar.get_parameters ( V0,nu0 );
        V0 *=1.0/nu0/10000.; // same flattening factor as for nu0
        nu0=1./10000.;*/
        Ar.set_parameters ( V0,nu0 );

        for ( int i=ilow; i<=ihi; i++ ) {
                Ar.set_parameters ( V0,nu0 );
                for ( int t=i; 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] );
        }

        if ( n_output>2 ) {
                Ar.set_parameters ( V0,nu0 );
                // Redo Ar for given points
                for ( int i=ilow; i<=ihi; i++ ) {
                        rgr.set_subvector ( 0,Data.get_col ( i ) );
                        Ar.bayes ( rgr );
                }

                // Create output vectors
                output[2] = mxCreateDoubleMatrix ( 3,1, mxREAL );
                // Convert the IT++ format to Matlab format for output
//              mat2mxArray ( Ar._epdf().mean(), output[2] );
                mat2mxArray ( Ar._epdf().mean(), output[2] );
        }

}