snipplrCallback({"title":"Mex Function for Mahalanobis Function","source":"
  1. #include "mex.h"<\/span><\/div><\/li>
  2. #include "math.h"<\/span><\/div><\/li>
  3. #include<armadillo> <\/span><\/div><\/li>
  4.  <\/div><\/li>
  5. using<\/span> namespace<\/span> arma;<\/div><\/li>
  6.  <\/div><\/li>
  7. void<\/span> importMatlab(<\/span>mat&<\/span> A, const<\/span> mxArray *<\/span>mxdata)<\/span>{<\/span><\/div><\/li>
  8.   access::<\/span>rw<\/span>(<\/span>A.mem<\/span>)<\/span>=<\/span>mxGetPr(<\/span>mxdata)<\/span>;<\/div><\/li>
  9.   access::<\/span>rw<\/span>(<\/span>A.n_rows<\/span>)<\/span>=<\/span>mxGetM(<\/span>mxdata)<\/span>;<\/div><\/li>
  10.   access::<\/span>rw<\/span>(<\/span>A.n_cols<\/span>)<\/span>=<\/span>mxGetN(<\/span>mxdata)<\/span>;<\/div><\/li>
  11.   access::<\/span>rw<\/span>(<\/span>A.n_elem<\/span>)<\/span>=<\/span>A.n_rows<\/span>*<\/span>A.n_cols<\/span>;<\/div><\/li>
  12. }<\/span>;<\/div><\/li>
  13.  <\/div><\/li>
  14. void<\/span> freeVar(<\/span>mat&<\/span> A, const<\/span> double<\/span> *<\/span>ptr)<\/span>{<\/span><\/div><\/li>
  15.   access::<\/span>rw<\/span>(<\/span>A.mem<\/span>)<\/span>=<\/span>ptr;<\/div><\/li>
  16.   access::<\/span>rw<\/span>(<\/span>A.n_rows<\/span>)<\/span>=<\/span>1<\/span>;<\/div><\/li>
  17.   access::<\/span>rw<\/span>(<\/span>A.n_cols<\/span>)<\/span>=<\/span>1<\/span>;<\/div><\/li>
  18.   access::<\/span>rw<\/span>(<\/span>A.n_elem<\/span>)<\/span>=<\/span>1<\/span>;<\/div><\/li>
  19. }<\/span>;<\/div><\/li>
  20.  <\/div><\/li>
  21.  <\/div><\/li>
  22. void<\/span> mexFunction(<\/span>int<\/span> nlhs, mxArray *<\/span>plhs[<\/span>]<\/span>, int<\/span> nrhs, const<\/span> mxArray *<\/span>prhs[<\/span>]<\/span>)<\/span><\/div><\/li>
  23. {<\/span><\/div><\/li>
  24.     if<\/span> (<\/span>nrhs !<\/span>=<\/span> 3<\/span>)<\/span><\/div><\/li>
  25.         mexErrMsgTxt(<\/span>"Incorrect number of input arguments"<\/span>)<\/span>;<\/div><\/li>
  26.     if<\/span> (<\/span>nlhs !<\/span>=<\/span> 1<\/span>)<\/span><\/div><\/li>
  27.         mexErrMsgTxt(<\/span>"Incorrect number of D arguments"<\/span>)<\/span>;<\/div><\/li>
  28.  <\/div><\/li>
  29.   mat X(<\/span>1<\/span>,1<\/span>)<\/span>;<\/div><\/li>
  30.   const<\/span> double<\/span>*<\/span> Xmem=<\/span>access::<\/span>rw<\/span>(<\/span>X.mem<\/span>)<\/span>;<\/div><\/li>
  31.   importMatlab(<\/span>X,prhs[<\/span>0<\/span>]<\/span>)<\/span>;<\/div><\/li>
  32.  <\/div><\/li>
  33.   mat Y(<\/span>1<\/span>,1<\/span>)<\/span>;<\/div><\/li>
  34.   const<\/span> double<\/span>*<\/span> Ymem=<\/span>access::<\/span>rw<\/span>(<\/span>Y.mem<\/span>)<\/span>;<\/div><\/li>
  35.   importMatlab(<\/span>Y,prhs[<\/span>1<\/span>]<\/span>)<\/span>;<\/div><\/li>
  36.  <\/div><\/li>
  37.   mat W(<\/span>1<\/span>,1<\/span>)<\/span>;<\/div><\/li>
  38.   const<\/span> double<\/span>*<\/span> Wmem=<\/span>access::<\/span>rw<\/span>(<\/span>W.mem<\/span>)<\/span>;<\/div><\/li>
  39.   importMatlab(<\/span>W,prhs[<\/span>2<\/span>]<\/span>)<\/span>;<\/div><\/li>
  40.  <\/div><\/li>
  41. \/\/    check if the input corresponds to what you are expecting<\/span><\/div><\/li>
  42.   if<\/span>(<\/span>(<\/span> X.n_cols<\/span> !<\/span>=<\/span>  Y.n_cols<\/span> )<\/span> || (<\/span>Y.n_cols<\/span> !<\/span>=<\/span> W.n_cols<\/span> )<\/span> )<\/span><\/div><\/li>
  43.     mexErrMsgTxt(<\/span>"Columns of X, Y, and W must be of equal length!"<\/span>)<\/span>;<\/div><\/li>
  44.    if<\/span>(<\/span> W.n_rows<\/span> !<\/span>=<\/span>  W.n_cols<\/span> )<\/span><\/div><\/li>
  45.     mexErrMsgTxt(<\/span>"W must be a square matrix!"<\/span>)<\/span>;<\/div><\/li>
  46.  <\/div><\/li>
  47.     plhs[<\/span>0<\/span>]<\/span> =<\/span> mxCreateDoubleMatrix(<\/span>X.n_rows<\/span>, Y.n_rows<\/span>, mxREAL)<\/span>;<\/div><\/li>
  48.  <\/div><\/li>
  49.     double<\/span> *<\/span>out =<\/span> mxGetPr(<\/span>plhs[<\/span>0<\/span>]<\/span>)<\/span>;<\/div><\/li>
  50.     mat diff;<\/div><\/li>
  51.     mat M;<\/div><\/li>
  52.     int<\/span> k=<\/span>0<\/span>;<\/div><\/li>
  53.     for<\/span>(<\/span>int<\/span> y=<\/span>0<\/span>;y<<\/span>Y.n_rows<\/span>;y++<\/span>)<\/span><\/div><\/li>
  54.       for<\/span>(<\/span>int<\/span> x=<\/span>0<\/span>;x<<\/span>X.n_rows<\/span>;x++<\/span>)<\/span>{<\/span><\/div><\/li>
  55. \/\/ (X(x,:)-Y(y,:))*W*((X(x,:)-Y(y,:))')<\/span><\/div><\/li>
  56.             diff=<\/span>X.row<\/span>(<\/span>x)<\/span>-<\/span>Y.row<\/span>(<\/span>y)<\/span>;<\/div><\/li>
  57.             M=<\/span>diff*<\/span>W*<\/span>trans(<\/span>diff)<\/span>;<\/div><\/li>
  58.             (<\/span>*<\/span>out++<\/span>)<\/span>=<\/span>M(<\/span>0<\/span>)<\/span>;<\/div><\/li>
  59.         }<\/span><\/div><\/li>
  60.  <\/div><\/li>
  61.  <\/div><\/li>
  62.     freeVar(<\/span>X,Xmem)<\/span>;<\/div><\/li>
  63.     freeVar(<\/span>Y,Ymem)<\/span>;<\/div><\/li>
  64.     freeVar(<\/span>W,Wmem)<\/span>;<\/div><\/li>
  65.     return<\/span>;<\/div><\/li>
  66. }<\/span><\/div><\/li><\/ol><\/pre>","link":"http:\/\/snipplr.com\/view\/22115\/mex-function-for-mahalanobis-function\/","id":"22115"});