idamnjanovic@51: /**************************************************************************
idamnjanovic@51:  *
idamnjanovic@51:  * File name: omp2mex.c
idamnjanovic@51:  *
idamnjanovic@51:  * Ron Rubinstein
idamnjanovic@51:  * Computer Science Department
idamnjanovic@51:  * Technion, Haifa 32000 Israel
idamnjanovic@51:  * ronrubin@cs
idamnjanovic@51:  *
idamnjanovic@51:  * Last Updated: 18.8.2009
idamnjanovic@51:  *
idamnjanovic@51:  *************************************************************************/
idamnjanovic@51: 
idamnjanovic@51: #include "ompcore.h"
idamnjanovic@51: #include "omputils.h"
idamnjanovic@51: #include "mexutils.h"
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* Input Arguments */
idamnjanovic@51: 
idamnjanovic@51: #define	IN_D	        prhs[0]
idamnjanovic@51: #define IN_X          prhs[1]
idamnjanovic@51: #define IN_DtX        prhs[2]
idamnjanovic@51: #define IN_XtX        prhs[3]
idamnjanovic@51: #define IN_G          prhs[4]
idamnjanovic@51: #define IN_EPS        prhs[5]
idamnjanovic@51: #define IN_SPARSE_G   prhs[6]
idamnjanovic@51: #define IN_MSGDELTA   prhs[7]
idamnjanovic@51: #define IN_MAXATOMS   prhs[8]
idamnjanovic@51: #define IN_PROFILE    prhs[9]
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* Output Arguments */
idamnjanovic@51: 
idamnjanovic@51: #define	GAMMA_OUT     plhs[0]
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /***************************************************************************************/
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[])
idamnjanovic@51: 
idamnjanovic@51: {
idamnjanovic@51:   double *D, *x, *DtX, *XtX, *G, eps, msgdelta;
idamnjanovic@51:   int gmode, maxatoms, profile;
idamnjanovic@51:   mwSize m, n, L;    /* D is n x m , X is n x L, DtX is m x L */
idamnjanovic@51: 
idamnjanovic@51:   
idamnjanovic@51:   /* check parameters */
idamnjanovic@51:   
idamnjanovic@51:   checkmatrix(IN_D, "OMP2", "D");
idamnjanovic@51:   checkmatrix(IN_X, "OMP2", "X");
idamnjanovic@51:   checkmatrix(IN_DtX, "OMP2", "DtX");
idamnjanovic@51:   checkmatrix(IN_XtX, "OMP2", "XtX");
idamnjanovic@51:   checkmatrix(IN_G, "OMP2", "G");
idamnjanovic@51:   
idamnjanovic@51:   checkscalar(IN_EPS, "OMP2", "EPSILON");
idamnjanovic@51:   checkscalar(IN_SPARSE_G, "OMP2", "sparse_g");
idamnjanovic@51:   checkscalar(IN_MSGDELTA, "OMP2", "msgdelta");
idamnjanovic@51:   checkscalar(IN_MAXATOMS, "OMP2", "maxatoms");
idamnjanovic@51:   checkscalar(IN_PROFILE, "OMP2", "profile");
idamnjanovic@51:   
idamnjanovic@51:   
idamnjanovic@51:   /* get parameters */
idamnjanovic@51:   
idamnjanovic@51:   x = D = DtX = XtX = G = 0;
idamnjanovic@51:   
idamnjanovic@51:   if (!mxIsEmpty(IN_D))
idamnjanovic@51:     D = mxGetPr(IN_D);
idamnjanovic@51:   
idamnjanovic@51:   if (!mxIsEmpty(IN_X))
idamnjanovic@51:     x = mxGetPr(IN_X);
idamnjanovic@51:   
idamnjanovic@51:   if (!mxIsEmpty(IN_DtX))
idamnjanovic@51:     DtX = mxGetPr(IN_DtX);
idamnjanovic@51:   
idamnjanovic@51:   if (!mxIsEmpty(IN_XtX))
idamnjanovic@51:     XtX = mxGetPr(IN_XtX);
idamnjanovic@51:   
idamnjanovic@51:   if (!mxIsEmpty(IN_G))
idamnjanovic@51:     G = mxGetPr(IN_G);
idamnjanovic@51:   
idamnjanovic@51:   eps = mxGetScalar(IN_EPS);
idamnjanovic@51:   if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) {
idamnjanovic@51:     gmode = SPARSE_GAMMA;
idamnjanovic@51:   }
idamnjanovic@51:   else {
idamnjanovic@51:     gmode = FULL_GAMMA;
idamnjanovic@51:   }
idamnjanovic@51:   msgdelta = mxGetScalar(IN_MSGDELTA);
idamnjanovic@51:   if (mxGetScalar(IN_MAXATOMS) < -1e-5) {
idamnjanovic@51:     maxatoms = -1;
idamnjanovic@51:   }
idamnjanovic@51:   else {
idamnjanovic@51:     maxatoms = (int)(mxGetScalar(IN_MAXATOMS)+1e-2);
idamnjanovic@51:   }
idamnjanovic@51:   profile = (int)(mxGetScalar(IN_PROFILE)+1e-2);
idamnjanovic@51:   
idamnjanovic@51:   
idamnjanovic@51:   /* check sizes */
idamnjanovic@51:   
idamnjanovic@51:   if (D && x) {
idamnjanovic@51:     n = mxGetM(IN_D);
idamnjanovic@51:     m = mxGetN(IN_D);
idamnjanovic@51:     L = mxGetN(IN_X);
idamnjanovic@51:     
idamnjanovic@51:     if (mxGetM(IN_X) != n) {
idamnjanovic@51:       mexErrMsgTxt("D and X have incompatible sizes.");
idamnjanovic@51:     }
idamnjanovic@51:     
idamnjanovic@51:     if (G) {
idamnjanovic@51:       if (mxGetN(IN_G)!=mxGetM(IN_G)) {
idamnjanovic@51:         mexErrMsgTxt("G must be a square matrix.");
idamnjanovic@51:       }
idamnjanovic@51:       if (mxGetN(IN_G) != m) {
idamnjanovic@51:         mexErrMsgTxt("D and G have incompatible sizes.");
idamnjanovic@51:       }
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   else if (DtX && XtX) {
idamnjanovic@51:     m = mxGetM(IN_DtX);
idamnjanovic@51:     L = mxGetN(IN_DtX);
idamnjanovic@51:     
idamnjanovic@51:     /* set n to an arbitrary value that is at least the max possible number of selected atoms */
idamnjanovic@51:     
idamnjanovic@51:     if (maxatoms>0) {
idamnjanovic@51:       n = maxatoms;
idamnjanovic@51:     }
idamnjanovic@51:     else {
idamnjanovic@51:       n = m;
idamnjanovic@51:     }
idamnjanovic@51:     
idamnjanovic@51:     if ( !(mxGetM(IN_XtX)==L && mxGetN(IN_XtX)==1) && !(mxGetM(IN_XtX)==1 && mxGetN(IN_XtX)==L) ) {
idamnjanovic@51:       mexErrMsgTxt("DtX and XtX have incompatible sizes.");
idamnjanovic@51:     }
idamnjanovic@51:     
idamnjanovic@51:     if (mxGetN(IN_G)!=mxGetM(IN_G)) {
idamnjanovic@51:       mexErrMsgTxt("G must be a square matrix.");
idamnjanovic@51:     }
idamnjanovic@51:     if (mxGetN(IN_G) != m) {
idamnjanovic@51:       mexErrMsgTxt("DtX and G have incompatible sizes.");
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   else {
idamnjanovic@51:     mexErrMsgTxt("Either D and X, or DtX and XtX, must be specified.");
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   
idamnjanovic@51:   /* Do OMP! */
idamnjanovic@51:   
idamnjanovic@51:   GAMMA_OUT = ompcore(D, x, DtX, XtX, G, n, m, L, maxatoms, eps, gmode, profile, msgdelta, 1);
idamnjanovic@51:   
idamnjanovic@51:   return;
idamnjanovic@51: }