finish adapters ; almost finished sources
authorBenjamin Auder <benjamin.auder@somewhere>
Thu, 15 Dec 2016 12:17:04 +0000 (13:17 +0100)
committerBenjamin Auder <benjamin.auder@somewhere>
Thu, 15 Dec 2016 12:17:04 +0000 (13:17 +0100)
src/adapters/a.EMGrank.c
src/adapters/a.constructionModelesLassoMLE.c
src/adapters/a.constructionModelesLassoRank.c
src/adapters/a.selectiontotale.c
src/sources/EMGrank.c

index 778ef67..f7fc192 100644 (file)
@@ -1,90 +1,72 @@
-#include "ioutils.h"
-#include "EMGrank.h"
-#include <mex.h>
-
-// nlhs, nrhs: resp. numbers of out and in parameters.
-// plhs: array of out parameters, each being a mxArray
-// plhs: array of in parameters (immutable), each being a mxArray
-//
-// MATLAB translates a call [A,B] = fun(C,D) into mexFunction(2,{A,B},2,{C,D}). 
-// Then mxArrayS are adapted to be passed to a regular C function, 
-// and the results are translated back to mxArrayS into plhs.
-void mexFunction(
-       int nlhs, 
-       mxArray* plhs[], 
-       int nrhs, 
-       const mxArray* prhs[])
-{
-       // Basic sanity checks
-       if (nrhs!=8) 
-               mexErrMsgIdAndTxt("select:EMGrank:nrhs","8 inputs required.");
-       if (nlhs!=2) 
-               mexErrMsgIdAndTxt("select:EMGrank:nlhs","3 outputs required.");
+#include <R.h>
+#include <Rdefines.h>
+#include "sources/EMGLLF.h"
 
+SEXP EMGLLF(
+       SEXP Pi_,
+       SEXP Rho_,
+       SEXP mini_,
+       SEXP maxi_,
+       SEXP X_,
+       SEXP Y_,
+       SEXP tau_,
+       SEXP rank_
+) {
        // Get matrices dimensions
-       const mwSize n = mxGetDimensions(prhs[4])[0];
-       const mwSize p = mxGetDimensions(prhs[4])[1];
-       const mwSize m = mxGetDimensions(prhs[1])[0];
-       const mwSize k = mxGetDimensions(prhs[1])[2];
+       SEXP dimX = getAttrib(X_, R_DimSymbol);
+       int n = INTEGER(dimX)[0];
+       int p = INTEGER(dimX)[1];
+       SEXP dimRho = getAttrib(Rho_, R_DimSymbol)
+       int m = INTEGER(dimRho)[0];
+       int k = INTEGER(dimRho)[2];
 
        ////////////
        // INPUTS //
        ////////////
 
-       // Pi
-       Real* Pi = mxGetPr(prhs[0]);
-
-       // Rho
-       const mwSize* dimRho = mxGetDimensions(prhs[1]);
-       Real* brRho = matlabToBrArray_real(mxGetPr(prhs[1]), dimRho, 3);
-       
-       // min number of iterations
-       Int mini = ((Int*)mxGetData(prhs[2]))[0];
-
-       // max number of iterations
-       Int maxi = ((Int*)mxGetData(prhs[3]))[0];
+       // get scalar parameters
+       int mini = INTEGER_VALUE(mini_);
+       int maxi = INTEGER_VALUE(maxi_);
+       double tau = NUMERIC_VALUE(tau_);
 
-       // X
-       const mwSize* dimX = mxGetDimensions(prhs[4]);
-       Real* brX = matlabToBrArray_real(mxGetPr(prhs[4]), dimX, 2);
-       
-       // Y
-       const mwSize* dimY = mxGetDimensions(prhs[5]);
-       Real* brY = matlabToBrArray_real(mxGetPr(prhs[5]), dimY, 2);
+       // Get pointers from SEXP arrays ; WARNING: by columns !
+       double* Pi = REAL(Pi_);
+       double* Rho = REAL(Rho_);
+       double* X = REAL(X_);
+       double* Y = REAL(Y_);
+       double* rank = REAL(rank_);
 
-       // tau
-       Real tau = mxGetScalar(prhs[6]);
-
-       // rank
-       Int* rank = (Int*)mxGetData(prhs[7]);
-       
        /////////////
        // OUTPUTS //
        /////////////
 
-       // phi
-       const mwSize dimPhi[] = {p, m, k};
-       plhs[0] = mxCreateNumericArray(3,dimPhi,mxDOUBLE_CLASS,mxREAL);
-       Real* phi = mxGetPr(plhs[0]);
-
-       // LLF
-       plhs[1] = mxCreateNumericMatrix(1,1,mxDOUBLE_CLASS,mxREAL);
-       Real* LLF = mxGetPr(plhs[1]);
+       SEXP phi, LLF;
+       PROTECT(dimPhi = allocVector(INTSXP, 3));
+       int* pDimPhi = INTEGER(dimPhi);
+       pDimPhi[0] = p; pDimPhi[1] = m; pDimPhi[2] = k;
+       PROTECT(phi = allocArray(REALSXP, dimPhi));
+       PROTECT(LLF = allocVector(REALSXP, 1));
+       double* pPhi=REAL(phi), pLLF=REAL(LLF);
 
        /////////////////////
        // Call to EMGrank //
        /////////////////////
 
-       EMGrank(Pi,brRho,mini,maxi,brX,brY,tau,rank,
-               phi,LLF,
+       EMGrank(Pi, Rho, mini, maxi, X, Y, tau, rank,
+               pPhi,pLLF,
                n,p,m,k);
 
-       free(brRho);
-       free(brX);
-       free(brY);
-       
-       //post-processing: convert by-rows outputs to MATLAB matrices
-       Real* mlPhi = brToMatlabArray_real(phi, dimPhi, 3);
-       copyArray(mlPhi, phi, dimPhi[0]*dimPhi[1]*dimPhi[2]);
-       free(mlPhi);
+       // Build list from OUT params and return it
+       SEXP listParams, listNames;
+       PROTECT(listParams = allocVector(VECSXP, 2));
+       char* lnames[2] = {"phi", "LLF"}; //lists labels
+       PROTECT(listNames = allocVector(STRSXP,2));
+       for (int i=0; i<2; i++)
+               SET_STRING_ELT(listNames,i,mkChar(lnames[i]));
+       setAttrib(listParams, R_NamesSymbol, listNames);
+       SET_ARRAY_ELT(listParams, 0, phi);
+       SET_VECTOR_ELT(listParams, 1, LLF);
+
+       UNPROTECT(5);
+       return listParams;
 }
index a331c67..7578664 100644 (file)
-#include "ioutils.h"
-#include "constructionModelesLassoMLE.h"
-#include <mex.h>
-
-#include <stdio.h>
-
-// nlhs, nrhs: resp. numbers of out and in parameters.
-// plhs: array of out parameters, each being a mxArray
-// plhs: array of in parameters (immutable), each being a mxArray
-//
-// MATLAB translates a call [A,B] = fun(C,D) into mexFunction(2,{A,B},2,{C,D}). 
-// Then mxArrayS are adapted to be passed to a regular C function, 
-// and the results are translated back to mxArrayS into plhs.
-void mexFunction(
-       int nlhs, 
-       mxArray* plhs[],
-       int nrhs, 
-       const mxArray* prhs[])
-{
-       // Basic sanity checks
-       if (nrhs!=14) 
-               mexErrMsgIdAndTxt("select:constructionModelesLassoMLE:nrhs","14 inputs required.");
-       if (nlhs!=4) 
-               mexErrMsgIdAndTxt("select:constructionModelesLassoMLE:nlhs","4 outputs required.");
-
+#include <R.h>
+#include <Rdefines.h>
+#include "sources/EMGLLF.h"
+
+SEXP EMGLLF(
+       SEXP phiInit_,
+       SEXP rhoInit_,
+       SEXP piInit_,
+       SEXP gamInit_,
+       SEXP mini_,
+       SEXP maxi_,
+       SEXP gamma_,
+       SEXP glambda_,
+       SEXP X_,
+       SEXP Y_,
+       SEXP seuil_,
+       SEXP tau_,
+       SEXP A1_,
+       SEXP A2_
+) {
        // Get matrices dimensions
-       const mwSize n = mxGetDimensions(prhs[8])[0];
-       const mwSize p = mxGetDimensions(prhs[0])[0];
-       const mwSize m = mxGetDimensions(prhs[0])[1];
-       const mwSize k = mxGetDimensions(prhs[0])[2];
-       const mwSize L = mxGetNumberOfElements(prhs[7]);
-       
+       int n = INTEGER(getAttrib(X_, R_DimSymbol))[0];
+       SEXP dim = getAttrib(phiInit_, R_DimSymbol)
+       int p = INTEGER(dim)[0];
+       int m = INTEGER(dim)[1];
+       int k = INTEGER(dim)[2];
+       int L = INTEGER(getAttrib(glambda_, R_LengthSymbol))[0];
+
        ////////////
        // INPUTS //
        ////////////
 
-       // phiInit
-       const mwSize* dimPhiInit = mxGetDimensions(prhs[0]);
-       Real* brPhiInit = matlabToBrArray_real(mxGetPr(prhs[0]), dimPhiInit, 3);
-       
-       // rhoInit
-       const mwSize* dimRhoInit = mxGetDimensions(prhs[1]);
-       Real* brRhoInit = matlabToBrArray_real(mxGetPr(prhs[1]), dimRhoInit, 3);
-
-       // piInit
-       Real* piInit = mxGetPr(prhs[2]);
-
-       // gamInit
-       const mwSize* dimGamInit = mxGetDimensions(prhs[3]);
-       Real* brGamInit = matlabToBrArray_real(mxGetPr(prhs[3]), dimGamInit, 2);
-
-       // min number of iterations
-       Int mini = ((Int*)mxGetData(prhs[4]))[0];
+       // get scalar parameters
+       int mini = INTEGER_VALUE(mini_);
+       int maxi = INTEGER_VALUE(maxi_);
+       double gamma = NUMERIC_VALUE(gamma_);
+       double seuil = NUMERIC_VALUE(seuil_);
+       double tau = NUMERIC_VALUE(tau_);
+
+       // Get pointers from SEXP arrays ; WARNING: by columns !
+       double* phiInit = REAL(phiInit_);
+       double* rhoInit = REAL(rhoInit_);
+       double* piInit = REAL(piInit_);
+       double* gamInit = REAL(gamInit_);
+       double* glambda = REAL(glambda_);
+       double* X = REAL(X_);
+       double* Y = REAL(Y_);
+       double* A1 = REAL(A1_);
+       double* A2 = REAL(A2_);
 
-       // max number of iterations
-       Int maxi = ((Int*)mxGetData(prhs[5]))[0];
-
-       // gamma
-       Real gamma = mxGetScalar(prhs[6]);
-
-       // glambda
-       Real* glambda = mxGetPr(prhs[7]);
-
-       // X
-       const mwSize* dimX = mxGetDimensions(prhs[8]);
-       Real* brX = matlabToBrArray_real(mxGetPr(prhs[8]), dimX, 2);
-
-       // Y
-       const mwSize* dimY = mxGetDimensions(prhs[9]);
-       Real* brY = matlabToBrArray_real(mxGetPr(prhs[9]), dimY, 2);
-       
-       //seuil
-       Real seuil = mxGetScalar(prhs[10]);
-
-       // tau
-       Real tau = mxGetScalar(prhs[11]);
-       
-       // A1
-       const mwSize* dimA = mxGetDimensions(prhs[12]);
-       Int* brA1 = matlabToBrArray_int(mxGetData(prhs[12]), dimA, 3);
-       
-       // A2
-       Int* brA2 = matlabToBrArray_int(mxGetData(prhs[13]), dimA, 3);
-       
        /////////////
        // OUTPUTS //
        /////////////
-       
-       // phi
-       const mwSize dimPhi[] = {dimPhiInit[0], dimPhiInit[1], dimPhiInit[2], L};
-       plhs[0] = mxCreateNumericArray(4,dimPhi,mxDOUBLE_CLASS,mxREAL);
-       Real* phi = mxGetPr(plhs[0]);
-       
-       // rho
-       const mwSize dimRho[] = {dimRhoInit[0], dimRhoInit[1], dimRhoInit[2], L};
-       plhs[1] = mxCreateNumericArray(4,dimRho,mxDOUBLE_CLASS,mxREAL);
-       Real* rho = mxGetPr(plhs[1]);
-
-       // pi
-       const mwSize dimPi[] = {k, L};
-       plhs[2] = mxCreateNumericMatrix(dimPi[0],dimPi[1],mxDOUBLE_CLASS,mxREAL);
-       Real* pi = mxGetPr(plhs[2]);
 
-       // lvraisemblance
-       const mwSize dimLvraisemblance[] = {L, 2};
-       plhs[3] = mxCreateNumericMatrix(L, 2, mxDOUBLE_CLASS,mxREAL);
-       Real* lvraisemblance = mxGetPr(plhs[3]);
+       SEXP phi, rho, pi, lvraisemblance, dimPhi, dimRho;
+       PROTECT(dimPhi = allocVector(INTSXP, 4));
+       int* pDimPhi = INTEGER(dimPhi);
+       pDimPhi[0] = p; pDimPhi[1] = m; pDimPhi[2] = k; pDimPhi[3] = L;
+       PROTECT(dimRho = allocVector(INTSXP, 4));
+       int* pDimRho = INTEGER(dimRho);
+       pDimRho[0] = m; pDimRho[1] = m; pDimRho[2] = k; pDimRho[3] = L;
+       PROTECT(phi = allocArray(REALSXP, dimPhi));
+       PROTECT(rho = allocArray(REALSXP, dimRho));
+       PROTECT(pi = allocMatrix(REALSXP, k, L));
+       PROTECT(lvraisemblance = allocMatrix(REALSXP, L, 2));
+       double* pPhi=REAL(phi), pRho=REAL(rho), pPi=REAL(pi), pLvraisemblance=REAL(lvraisemblance);
 
        /////////////////////////////////////////
        // Call to constructionModelesLassoMLE //
        /////////////////////////////////////////
 
        constructionModelesLassoMLE(
-               brPhiInit,brRhoInit,piInit,brGamInit,mini,maxi,gamma,glambda,brX,brY,seuil,tau,brA1,brA2,
-               phi,rho,pi,lvraisemblance,
+               phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2,
+               pPhi,pRho,pPi,pLvraisemblance,
                n,p,m,k,L);
-       
-       free(brPhiInit);
-       free(brRhoInit);
-       free(brGamInit);
-       free(brX);
-       free(brY);
-       free(brA1);
-       free(brA2);
-       
-       //post-processing: convert by-rows outputs to MATLAB matrices
-       Real* mlPhi = brToMatlabArray_real(phi, dimPhi, 4);
-       copyArray(mlPhi, phi, dimPhi[0]*dimPhi[1]*dimPhi[2]*dimPhi[3]);
-       free(mlPhi);
-       
-       Real* mlRho = brToMatlabArray_real(rho, dimRho, 4);
-       copyArray(mlRho, rho, dimRho[0]*dimRho[1]*dimRho[2]*dimRho[3]);
-       free(mlRho);
-
-       Real* mlPi = brToMatlabArray_real(pi, dimPi, 2);
-       copyArray(mlPi, pi, dimPi[0]*dimPi[1]);
-       free(mlPi);
 
-       Real* mlLvraisemblance = brToMatlabArray_real(lvraisemblance, dimLvraisemblance, 2);
-       copyArray(mlLvraisemblance, lvraisemblance, dimLvraisemblance[0]*dimLvraisemblance[1]);
-       free(mlLvraisemblance);
+       // Build list from OUT params and return it
+       SEXP listParams, listNames;
+       PROTECT(listParams = allocVector(VECSXP, 4));
+       char* lnames[4] = {"phi", "rho", "pi", "lvraisemblance}; //lists labels
+       PROTECT(listNames = allocVector(STRSXP,4));
+       for (int i=0; i<4; i++)
+               SET_STRING_ELT(listNames,i,mkChar(lnames[i]));
+       setAttrib(listParams, R_NamesSymbol, listNames);
+       SET_ARRAY_ELT(listParams, 0, phi);
+       SET_ARRAY_ELT(listParams, 1, rho);
+       SET_MATRIX_ELT(listParams, 2, pi);
+       SET_VECTOR_ELT(listParams, 3, lvraisemblance);
+
+       UNPROTECT(8);
+       return listParams;
 }
index 96908c2..359290c 100644 (file)
-#include "ioutils.h"
-#include "constructionModelesLassoRank.h"
-#include <mex.h>
+#include <R.h>
+#include <Rdefines.h>
+#include "sources/EMGLLF.h"
 
-#include <stdio.h>
-
-// nlhs, nrhs: resp. numbers of out and in parameters.
-// plhs: array of out parameters, each being a mxArray
-// plhs: array of in parameters (immutable), each being a mxArray
-//
-// MATLAB translates a call [A,B] = fun(C,D) into mexFunction(2,{A,B},2,{C,D}). 
-// Then mxArrayS are adapted to be passed to a regular C function, 
-// and the results are translated back to mxArrayS into plhs.
-void mexFunction(int nlhs, mxArray* plhs[], 
-                 int nrhs, const mxArray* prhs[])
-{
-       // Basic sanity checks
-       if (nrhs!=10) 
-               mexErrMsgIdAndTxt("select:constructionModelesLassoRank:nrhs","10 inputs required.");
-       if (nlhs!=2) 
-               mexErrMsgIdAndTxt("select:constructionModelesLassoRank:nlhs","3 outputs required.");
-
-       // Get matrices dimensions, to be given to main routine above
-       const mwSize n = mxGetDimensions(prhs[4])[0];
-       const mwSize p = mxGetDimensions(prhs[4])[1];
-       const mwSize m = mxGetDimensions(prhs[1])[0];
-       const mwSize k = mxGetDimensions(prhs[1])[2];
-       const mwSize L = mxGetDimensions(prhs[7])[1];
+SEXP EMGLLF(
+       SEXP Pi_,
+       SEXP Rho_,
+       SEXP mini_,
+       SEXP maxi_,
+       SEXP X_,
+       SEXP Y_,
+       SEXP tau_,
+       SEXP A1_,
+       SEXP rangmin_,
+       SEXP rangmax
+) {
+       // Get matrices dimensions
+       SEXP dimX = getAttrib(X_, R_DimSymbol);
+       int n = INTEGER(dimX)[0];
+       int p = INTEGER(dimX)[1];
+       SEXP dimRho = getAttrib(Rho_, R_DimSymbol)
+       int m = INTEGER(dimRho)[0];
+       int k = INTEGER(dimRho)[2];
+       int L = INTEGER(getAttrib(A1_, R_DimSymbol))[1];
 
        ////////////
        // INPUTS //
        ////////////
 
-       // pi
-       const mwSize* dimPi = mxGetDimensions(prhs[0]);
-       Real* brPi = matlabToBrArray_real(mxGetPr(prhs[0]), dimPi, 2);
-
-       // rho
-       const mwSize* dimRho = mxGetDimensions(prhs[1]);
-       Real* brRho = matlabToBrArray_real(mxGetPr(prhs[1]), dimRho, 4);
-
-       // min number of iterations
-       Int mini = ((Int*)mxGetData(prhs[2]))[0];
-
-       // max number of iterations
-       Int maxi = ((Int*)mxGetData(prhs[3]))[0];
+       // get scalar parameters
+       int mini = INTEGER_VALUE(mini_);
+       int maxi = INTEGER_VALUE(maxi_);
+       double tau = NUMERIC_VALUE(tau_);
+       double rangmin = NUMERIC_VALUE(rangmin_);
+       double rangmax = NUMERIC_VALUE(rangmax_);
 
-       // X
-       const mwSize* dimX = mxGetDimensions(prhs[4]);
-       Real* brX = matlabToBrArray_real(mxGetPr(prhs[4]), dimX, 2);
+       // Get pointers from SEXP arrays ; WARNING: by columns !
+       double* Pi = REAL(Pi_);
+       double* Rho = REAL(Rho_);
+       double* X = REAL(X_);
+       double* Y = REAL(Y_);
+       double* A1 = REAL(A1_);
 
-       // Y
-       const mwSize* dimY = mxGetDimensions(prhs[5]);
-       Real* brY = matlabToBrArray_real(mxGetPr(prhs[5]), dimY, 2);
-       
-       // tau
-       Real tau = mxGetScalar(prhs[6]);
-       
-       // A1
-       const mwSize* dimA = mxGetDimensions(prhs[7]);
-       Int* brA1 = matlabToBrArray_int(mxGetData(prhs[7]), dimA, 2);
-       
-    //rangmin
-    Int rangmin = ((Int*)mxGetData(prhs[8]))[0];
-    
-    //rangmax
-    Int rangmax = ((Int*)mxGetData(prhs[9]))[0];
-       
        /////////////
        // OUTPUTS //
        /////////////
-       
-       // phi
-       mwSize Size = pow(rangmax-rangmin+1,k);
-       const mwSize dimPhi[] = {p, m, k, L*Size};
-       plhs[0] = mxCreateNumericArray(4,dimPhi,mxDOUBLE_CLASS,mxREAL);
-       Real* phi = mxGetPr(plhs[0]);
 
-       // lvraisemblance
-       const mwSize dimLvraisemblance[] = {L*Size, 2};
-       plhs[1] = mxCreateNumericMatrix(dimLvraisemblance[0],dimLvraisemblance[1],mxDOUBLE_CLASS,mxREAL);
-       Real* lvraisemblance = mxGetPr(plhs[1]);
+       int Size = pow(rangmax-rangmin+1,k);
+       SEXP phi, lvraisemblance, dimPhi;
+       PROTECT(dimPhi = allocVector(INTSXP, 4));
+       int* pDimPhi = INTEGER(dimPhi);
+       pDimPhi[0] = p; pDimPhi[1] = m; pDimPhi[2] = k; pDimPhi[3] = L*Size;
+       PROTECT(phi = allocArray(REALSXP, dimPhi));
+       PROTECT(lvraisemblance = allocMatrix(REALSXP, L*Size, 2));
+       double* pPhi=REAL(phi), pLvraisemblance=REAL(lvraisemblance);
 
-       
        //////////////////////////////////////////
        // Call to constructionModelesLassoRank //
        //////////////////////////////////////////
 
        constructionModelesLassoRank(
-        brPi,brRho,mini,maxi,brX,brY,tau,brA1,rangmin,rangmax,
-        phi,lvraisemblance,
-        n,p,m,k,L);
-       
-       free(brPi);
-       free(brRho);
-       free(brX);
-       free(brY);
-       free(brA1);
+               Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax,
+               pPhi,pLvraisemblance,
+               n,p,m,k,L);
 
-       //post-processing: convert by-rows outputs to MATLAB matrices
-       Real* mlPhi = brToMatlabArray_real(phi, dimPhi, 4);
-       copyArray(mlPhi, phi, dimPhi[0]*dimPhi[1]*dimPhi[2]*dimPhi[3]);
-       free(mlPhi);
-       
-       Real* mlLvraisemblance = brToMatlabArray_real(lvraisemblance, dimLvraisemblance, 2);
-       copyArray(mlLvraisemblance, lvraisemblance, dimLvraisemblance[0]*dimLvraisemblance[1]);
-       free(mlLvraisemblance);
+       // Build list from OUT params and return it
+       SEXP listParams, listNames;
+       PROTECT(listParams = allocVector(VECSXP, 2));
+       char* lnames[2] = {"phi", "lvraisemblance"}; //lists labels
+       PROTECT(listNames = allocVector(STRSXP,2));
+       for (int i=0; i<2; i++)
+               SET_STRING_ELT(listNames,i,mkChar(lnames[i]));
+       setAttrib(listParams, R_NamesSymbol, listNames);
+       SET_ARRAY_ELT(listParams, 0, phi);
+       SET_VECTOR_ELT(listParams, 1, lvraisemblance);
 
+       UNPROTECT(5);
+       return listParams;
 }
index b4b8f4c..1124d71 100644 (file)
-#include "ioutils.h"
-#include "selectiontotale.h"
-#include "EMGLLF.h"
-#include <mex.h>
-
-// nlhs, nrhs: resp. numbers of out and in parameters.
-// plhs: array of out parameters, each being a mxArray
-// plhs: array of in parameters (immutable), each being a mxArray
-//
-// MATLAB translates a call [A,B] = fun(C,D) into mexFunction(2,{A,B},2,{C,D}). 
-// Then mxArrayS are adapted to be passed to a regular C function, 
-// and the results are translated back to mxArrayS into plhs.
-void mexFunction(int nlhs, mxArray* plhs[], 
-                 int nrhs, const mxArray* prhs[])
-{
-       // Basic sanity checks
-       if (nrhs!=12) 
-               mexErrMsgIdAndTxt("select:selectiontotale:nrhs","12 inputs required.");
-       if (nlhs!=4) 
-               mexErrMsgIdAndTxt("select:selectiontotale:nlhs","4 outputs required.");
-
-       // Get matrices dimensions, to be given to main routine above
-       const mwSize n = mxGetDimensions(prhs[8])[0];
-       const mwSize p = mxGetDimensions(prhs[8])[1];
-       const mwSize m = mxGetDimensions(prhs[1])[0];
-       const mwSize k = mxGetDimensions(prhs[1])[2];
-       const mwSize L = mxGetNumberOfElements(prhs[7]);
+#include <R.h>
+#include <Rdefines.h>
+#include "sources/EMGLLF.h"
+
+SEXP EMGLLF(
+       SEXP phiInit_,
+       SEXP rhoInit_,
+       SEXP piInit_,
+       SEXP gamInit_,
+       SEXP mini_,
+       SEXP maxi_,
+       SEXP gamma_,
+       SEXP glambda_,
+       SEXP X_,
+       SEXP Y_,
+       SEXP seuil_,
+       SEXP tau_
+) {
+       // Get matrices dimensions
+       SEXP dimX = getAttrib(X_, R_DimSymbol);
+       int n = INTEGER(dimX)[0];
+       int p = INTEGER(dimX)[1];
+       SEXP dimRho = getAttrib(rhoInit_, R_DimSymbol)
+       int m = INTEGER(dimRho)[0];
+       int k = INTEGER(dimRho)[2];
+       int L = INTEGER(getAttrib(glambda_, R_LengthSymbol))[0];
 
        ////////////
        // INPUTS //
        ////////////
 
-       // phiInit
-       const mwSize* dimPhiInit = mxGetDimensions(prhs[0]);
-       Real* brPhiInit = matlabToBrArray_real(mxGetPr(prhs[0]), dimPhiInit, 3);
-       
-       // rhoInit
-       const mwSize* dimRhoInit = mxGetDimensions(prhs[1]);
-       Real* brRhoInit = matlabToBrArray_real(mxGetPr(prhs[1]), dimRhoInit, 3);
-
-       // piInit
-       Real* piInit = mxGetPr(prhs[2]);
-
-       // gamInit
-       const mwSize* dimGamInit = mxGetDimensions(prhs[3]);
-       Real* brGamInit = matlabToBrArray_real(mxGetPr(prhs[3]), dimGamInit, 2);
-
-       // min number of iterations
-       Int mini = ((Int*)mxGetData(prhs[4]))[0];
-
-       // max number of iterations
-       Int maxi = ((Int*)mxGetData(prhs[5]))[0];
-
-       // gamma
-       Real gamma = mxGetScalar(prhs[6]);
-
-       // glambda
-       Real* glambda = mxGetPr(prhs[7]);
-
-       // X
-       const mwSize* dimX = mxGetDimensions(prhs[8]);
-       Real* brX = matlabToBrArray_real(mxGetPr(prhs[8]), dimX, 2);
-       
-       // Y
-       const mwSize* dimY = mxGetDimensions(prhs[9]);
-       Real* brY = matlabToBrArray_real(mxGetPr(prhs[9]), dimY, 2);
-       
-       //seuil
-       Real seuil = mxGetScalar(prhs[10]);
-
-       // tau
-       Real tau = mxGetScalar(prhs[11]);
+       // get scalar parameters
+       int mini = INTEGER_VALUE(mini_);
+       int maxi = INTEGER_VALUE(maxi_);
+       double gamma = NUMERIC_VALUE(gamma_);
+       double seuil = NUMERIC_VALUE(seuil_);
+       double tau = NUMERIC_VALUE(tau_);
+
+       // Get pointers from SEXP arrays ; WARNING: by columns !
+       double* piInit = REAL(phiInit_);
+       double* rhoInit = REAL(rhoInit_);
+       double* piInit = REAL(piInit_);
+       double* gamInit = REAL(gamInit_);
+       double* glambda = REAL(glambda_);
+       double* X = REAL(X_);
+       double* Y = REAL(Y_);
 
        /////////////
        // OUTPUTS //
        /////////////
-       
-       // A1
-       mwSize dimA[] = {p,m+1,L};
-       plhs[0] = mxCreateNumericArray(3,dimA,mxGetClassID(prhs[4]),mxREAL);
-       Int* A1 = (Int*)mxGetData(plhs[0]);
-       
-       // A2
-       plhs[1] = mxCreateNumericArray(3,dimA,mxGetClassID(prhs[4]),mxREAL);
-       Int* A2 = (Int*)mxGetData(plhs[1]);
-       
-       // rho
-       const mwSize dimRho[] = {dimRhoInit[0], dimRhoInit[1], dimRhoInit[2], L};
-       plhs[2] = mxCreateNumericArray(4,dimRho,mxDOUBLE_CLASS,mxREAL);
-       Real* Rho = mxGetPr(plhs[2]);
 
-       // pi
-       const mwSize dimPi[] = {k, L};
-       plhs[3] = mxCreateNumericMatrix(dimPi[0],dimPi[1],mxDOUBLE_CLASS,mxREAL);
-       double* Pi = mxGetPr(plhs[3]);
+       int Size = pow(rangmax-rangmin+1,k);
+       SEXP A1, A2, rho, pi, dimA;
+       PROTECT(dimA = allocVector(INTSXP, 3));
+       int* pDimA = INTEGER(dimA);
+       pDimA[0] = p; pDimA[1] = m+1; pDimA[2] = L;
+       PROTECT(A1 = allocArray(REALSXP, dimA));
+       PROTECT(A2 = allocArray(REALSXP, dimA));
+       PROTECT(rho = allocArray(REALSXP, dimRho);
+       PROTECT(pi = allocMatrix(REALSXP, k, L));
+       double* pA1=REAL(A1), pA2=REAL(A2), pRho=REAL(rho), pPi=REAL(pi);
 
        /////////////////////////////
        // Call to selectiontotale //
        /////////////////////////////
 
-       selectiontotale(brPhiInit,brRhoInit,piInit,brGamInit,mini,maxi,gamma,glambda,brX,brY,seuil,tau,
-               A1,A2,Rho,Pi,
+       selectiontotale(
+               phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,
+               pA1,pA2,pRho,pPi,
                n,p,m,k,L);
-       
-       free(brPhiInit);
-       free(brRhoInit);
-       free(brGamInit);
-       free(brX);
-       free(brY);
-       
-       //post-processing: convert by-rows outputs to MATLAB matrices
-       Int* mlA1 = brToMatlabArray_int(A1,dimA,3);
-       copyArray(mlA1,A1,dimA[0]*dimA[1]*dimA[2]);
-       free(mlA1);
-       
-       Int* mlA2 = brToMatlabArray_int(A2,dimA,3);
-       copyArray(mlA2,A2,dimA[0]*dimA[1]*dimA[2]);
-       free(mlA2);
-       
-       Real* mlRho = brToMatlabArray_real(Rho, dimRho, 4);
-       copyArray(mlRho, Rho, dimRho[0]*dimRho[1]*dimRho[2]*dimRho[3]);
-       free(mlRho);
 
-       Real* mlPi = brToMatlabArray_real(Pi, dimPi, 2);
-       copyArray(mlPi, Pi, dimPi[0]*dimPi[1]);
-       free(mlPi);
+       // Build list from OUT params and return it
+       SEXP listParams, listNames;
+       PROTECT(listParams = allocVector(VECSXP, 4));
+       char* lnames[4] = { "A1", "A2", "rho", "pi" }; //lists labels
+       PROTECT(listNames = allocVector(STRSXP, 4));
+       for (int i=0; i<4; i++)
+               SET_STRING_ELT(listNames,i,mkChar(lnames[i]));
+       setAttrib(listParams, R_NamesSymbol, listNames);
+       SET_ARRAY_ELT(listParams, 0, A1);
+       SET_ARRAY_ELT(listParams, 1, A2);
+       SET_ARRAY_ELT(listParams, 2, rho);
+       SET_MATRIX_ELT(listParams, 3, pi);
+
+       UNPROTECT(7);
+       return listParams;
 }
index 721c0c2..2e85ce6 100644 (file)
@@ -2,30 +2,29 @@
 #include <gsl/gsl_linalg.h>
 
 // Compute pseudo-inverse of a square matrix
-static Real* pinv(const Real* matrix, mwSize dim)
+static double* pinv(const double* matrix, int dim)
 {
        gsl_matrix* U = gsl_matrix_alloc(dim,dim);
        gsl_matrix* V = gsl_matrix_alloc(dim,dim);
        gsl_vector* S = gsl_vector_alloc(dim);
        gsl_vector* work = gsl_vector_alloc(dim);
-       Real EPS = 1e-10; //threshold for singular value "== 0"
+       double EPS = 1e-10; //threshold for singular value "== 0"
        
        //copy matrix into U
-       for (mwSize i=0; i<dim*dim; i++)
-               U->data[i] = matrix[i];
-       
+       copyArray(matrix, U->data, dim*dim);
+
        //U,S,V = SVD of matrix
        gsl_linalg_SV_decomp(U, V, S, work);
        gsl_vector_free(work);
-       
+
        // Obtain pseudo-inverse by V*S^{-1}*t(U)
-       Real* inverse = (Real*)malloc(dim*dim*sizeof(Real));
-       for (mwSize i=0; i<dim; i++)
+       double* inverse = (double*)malloc(dim*dim*sizeof(double));
+       for (int i=0; i<dim; i++)
        {
-               for (mwSize ii=0; ii<dim; ii++)
+               for (int ii=0; ii<dim; ii++)
                {
-                       Real dotProduct = 0.0;
-                       for (mwSize j=0; j<dim; j++)
+                       double dotProduct = 0.0;
+                       for (int j=0; j<dim; j++)
                                dotProduct += V->data[i*dim+j] * (S->data[j] > EPS ? 1.0/S->data[j] : 0.0) * U->data[ii*dim+j];
                        inverse[i*dim+ii] = dotProduct;
                }
@@ -40,39 +39,39 @@ static Real* pinv(const Real* matrix, mwSize dim)
 // TODO: comment EMGrank purpose
 void EMGrank(
        // IN parameters
-       const Real* Pi,          // parametre de proportion
-       const Real* Rho,         // parametre initial de variance renormalisé
-       Int mini,            // nombre minimal d'itérations dans l'algorithme EM        
-       Int maxi,            // nombre maximal d'itérations dans l'algorithme EM
-       const Real* X,           // régresseurs
-       const Real* Y,           // réponse
-       Real tau,          // seuil pour accepter la convergence
-       const Int* rank,           // vecteur des rangs possibles
+       const double* Pi, // parametre de proportion
+       const double* Rho, // parametre initial de variance renormalisé
+       int mini, // nombre minimal d'itérations dans l'algorithme EM
+       int maxi, // nombre maximal d'itérations dans l'algorithme EM
+       const double* X, // régresseurs
+       const double* Y, // réponse
+       double tau, // seuil pour accepter la convergence
+       const int* rank, // vecteur des rangs possibles
        // OUT parameters
-       Real* phi,         // parametre de moyenne renormalisé, calculé par l'EM
-       Real* LLF,         // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres
+       double* phi, // parametre de moyenne renormalisé, calculé par l'EM
+       double* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres
        // additional size parameters
-       mwSize n,               // taille de l'echantillon                
-       mwSize p,               // nombre de covariables
-       mwSize m,               // taille de Y (multivarié)
-       mwSize k)              // nombre de composantes
+       int n, // taille de l'echantillon
+       int p, // nombre de covariables
+       int m, // taille de Y (multivarié)
+       int k) // nombre de composantes
 {
        // Allocations, initializations
-       Real* Phi = (Real*)calloc(p*m*k,sizeof(Real));
-       Real* hatBetaR = (Real*)malloc(p*m*sizeof(Real));
+       double* Phi = (double*)calloc(p*m*k,sizeof(double));
+       double* hatBetaR = (double*)malloc(p*m*sizeof(double));
        int signum;
-       Real invN = 1.0/n;
+       double invN = 1.0/n;
        int deltaPhiBufferSize = 20;
-       Real* deltaPhi = (Real*)malloc(deltaPhiBufferSize*sizeof(Real));
-       mwSize ite = 0;
-       Real sumDeltaPhi = 0.0;
-       Real* YiRhoR = (Real*)malloc(m*sizeof(Real));
-       Real* XiPhiR = (Real*)malloc(m*sizeof(Real));
-       Real* Xr = (Real*)malloc(n*p*sizeof(Real));
-       Real* Yr = (Real*)malloc(n*m*sizeof(Real));
-       Real* tXrXr = (Real*)malloc(p*p*sizeof(Real));
-       Real* tXrYr = (Real*)malloc(p*m*sizeof(Real));
-       gsl_matrix* matrixM = gsl_matrix_alloc(p, m);        
+       double* deltaPhi = (double*)malloc(deltaPhiBufferSize*sizeof(double));
+       int ite = 0;
+       double sumDeltaPhi = 0.0;
+       double* YiRhoR = (double*)malloc(m*sizeof(double));
+       double* XiPhiR = (double*)malloc(m*sizeof(double));
+       double* Xr = (double*)malloc(n*p*sizeof(double));
+       double* Yr = (double*)malloc(n*m*sizeof(double));
+       double* tXrXr = (double*)malloc(p*p*sizeof(double));
+       double* tXrYr = (double*)malloc(p*m*sizeof(double));
+       gsl_matrix* matrixM = gsl_matrix_alloc(p, m);
        gsl_matrix* matrixE = gsl_matrix_alloc(m, m);
        gsl_permutation* permutation = gsl_permutation_alloc(m);
        gsl_matrix* V = gsl_matrix_alloc(m,m);
@@ -80,158 +79,158 @@ void EMGrank(
        gsl_vector* work = gsl_vector_alloc(m);
 
        //Initialize class memberships (all elements in class 0; TODO: randomize ?)
-       Int* Z = (Int*)calloc(n, sizeof(Int));
+       int* Z = (int*)calloc(n, sizeof(int));
 
        //Initialize phi to zero, because some M loops might exit before phi affectation
-       for (mwSize i=0; i<p*m*k; i++)
+       for (int i=0; i<p*m*k; i++)
                phi[i] = 0.0;
 
        while (ite<mini || (ite<maxi && sumDeltaPhi>tau))
-       {               
+       {
                /////////////
                // Etape M //
                /////////////
                
                //M step: Mise à jour de Beta (et donc phi)
-               for (mwSize r=0; r<k; r++)
+               for (int r=0; r<k; r++)
                {
                        //Compute Xr = X(Z==r,:) and Yr = Y(Z==r,:)
-                       mwSize cardClustR=0;
-                       for (mwSize i=0; i<n; i++)
+                       int cardClustR=0;
+                       for (int i=0; i<n; i++)
                        {
                                if (Z[i] == r)
                                {
-                                       for (mwSize j=0; j<p; j++)
-                                               Xr[cardClustR*p+j] = X[i*p+j];
-                                       for (mwSize j=0; j<m; j++)
-                                               Yr[cardClustR*m+j] = Y[i*m+j];
+                                       for (int j=0; j<p; j++)
+                                               Xr[mi(cardClustR,j,n,p)] = X[mi(i,j,n,p)];
+                                       for (int j=0; j<m; j++)
+                                               Yr[mi(cardClustR,j,n,m)] = Y[mi(i,j,n,m)];
                                        cardClustR++;
                                }
                        }
-                       if (cardClustR == 0) 
+                       if (cardClustR == 0)
                                continue;
 
                        //Compute tXrXr = t(Xr) * Xr
-                       for (mwSize j=0; j<p; j++)
+                       for (int j=0; j<p; j++)
                        {
-                               for (mwSize jj=0; jj<p; jj++)
+                               for (int jj=0; jj<p; jj++)
                                {
-                                       Real dotProduct = 0.0;
-                                       for (mwSize u=0; u<cardClustR; u++)
-                                               dotProduct += Xr[u*p+j] * Xr[u*p+jj];
-                                       tXrXr[j*p+jj] = dotProduct;
+                                       double dotProduct = 0.0;
+                                       for (int u=0; u<cardClustR; u++)
+                                               dotProduct += Xr[mi(u,j,n,p)] * Xr[mi(u,jj,n,p)];
+                                       tXrXr[mi(j,jj,p,p)] = dotProduct;
                                }
                        }
 
                        //Get pseudo inverse = (t(Xr)*Xr)^{-1}
-                       Real* invTXrXr = pinv(tXrXr, p);
+                       double* invTXrXr = pinv(tXrXr, p);
                        
                        // Compute tXrYr = t(Xr) * Yr
-                       for (mwSize j=0; j<p; j++)
+                       for (int j=0; j<p; j++)
                        {
-                               for (mwSize jj=0; jj<m; jj++)
+                               for (int jj=0; jj<m; jj++)
                                {
-                                       Real dotProduct = 0.0;
-                                       for (mwSize u=0; u<cardClustR; u++)
-                                               dotProduct += Xr[u*p+j] * Yr[u*m+jj];
-                                       tXrYr[j*m+jj] = dotProduct;
+                                       double dotProduct = 0.0;
+                                       for (int u=0; u<cardClustR; u++)
+                                               dotProduct += Xr[mi(u,j,n,p)] * Yr[mi(u,j,n,m)];
+                                       tXrYr[mi(j,jj,p,m)] = dotProduct;
                                }
                        }
 
                        //Fill matrixM with inverse * tXrYr = (t(Xr)*Xr)^{-1} * t(Xr) * Yr
-                       for (mwSize j=0; j<p; j++)
+                       for (int j=0; j<p; j++)
                        {
-                               for (mwSize jj=0; jj<m; jj++)
+                               for (int jj=0; jj<m; jj++)
                                {
-                                       Real dotProduct = 0.0;
-                                       for (mwSize u=0; u<p; u++)
-                                               dotProduct += invTXrXr[j*p+u] * tXrYr[u*m+jj];
-                                       matrixM->data[j*m+jj] = dotProduct;     
+                                       double dotProduct = 0.0;
+                                       for (int u=0; u<p; u++)
+                                               dotProduct += invTXrXr[mi(j,u,p,p)] * tXrYr[mi(u,jj,p,m)];
+                                       matrixM->data[j*m+jj] = dotProduct;
                                }
                        }
                        free(invTXrXr);
-                       
+
                        //U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr
                        gsl_linalg_SV_decomp(matrixM, V, S, work);
-                       
-                       //Set m-rank(r) singular values to zero, and recompose 
+
+                       //Set m-rank(r) singular values to zero, and recompose
                        //best rank(r) approximation of the initial product
-                       for (mwSize j=rank[r]; j<m; j++)
+                       for (int j=rank[r]; j<m; j++)
                                S->data[j] = 0.0;
                        
                        //[intermediate step] Compute hatBetaR = U * S * t(V)
-                       Real* U = matrixM->data;
-                       for (mwSize j=0; j<p; j++)
+                       double* U = matrixM->data;
+                       for (int j=0; j<p; j++)
                        {
-                               for (mwSize jj=0; jj<m; jj++)
+                               for (int jj=0; jj<m; jj++)
                                {
-                                       Real dotProduct = 0.0;
-                                       for (mwSize u=0; u<m; u++) 
+                                       double dotProduct = 0.0;
+                                       for (int u=0; u<m; u++)
                                                dotProduct += U[j*m+u] * S->data[u] * V->data[jj*m+u];
-                                       hatBetaR[j*m+jj] = dotProduct;
+                                       hatBetaR[mi(j,jj,p,m)] = dotProduct;
                                }
                        }
+
                        //Compute phi(:,:,r) = hatBetaR * Rho(:,:,r)
-                       for (mwSize j=0; j<p; j++)
+                       for (int j=0; j<p; j++)
                        {
-                               for (mwSize jj=0; jj<m; jj++)
+                               for (int jj=0; jj<m; jj++)
                                {
-                                       Real dotProduct=0.0;
-                                       for (mwSize u=0; u<m; u++)
-                                               dotProduct += hatBetaR[j*m+u] * Rho[u*m*k+jj*k+r];
-                                       phi[j*m*k+jj*k+r] = dotProduct;
+                                       double dotProduct=0.0;
+                                       for (int u=0; u<m; u++)
+                                               dotProduct += hatBetaR[mi(j,u,p,m)] * Rho[ai(u,jj,r,m,m,k)];
+                                       phi[ai(j,jj,r,p,m,k)] = dotProduct;
                                }
-                   }
+               }
                }
 
                /////////////
                // Etape E //
                /////////////
                
-               Real sumLogLLF2 = 0.0;
-               for (mwSize i=0; i<n; i++)
+               double sumLogLLF2 = 0.0;
+               for (int i=0; i<n; i++)
                {
-                       Real sumLLF1 = 0.0;
-                       Real maxLogGamIR = -INFINITY;
-                       for (mwSize r=0; r<k; r++)
+                       double sumLLF1 = 0.0;
+                       double maxLogGamIR = -INFINITY;
+                       for (int r=0; r<k; r++)
                        {
                                //Compute
                                //Gam(i,r) = Pi(r) * det(Rho(:,:,r)) * exp( -1/2 * (Y(i,:)*Rho(:,:,r) - X(i,:)...
-                               //    *phi(:,:,r)) * transpose( Y(i,:)*Rho(:,:,r) - X(i,:)*phi(:,:,r) ) );
+                               //*phi(:,:,r)) * transpose( Y(i,:)*Rho(:,:,r) - X(i,:)*phi(:,:,r) ) );
                                //split in several sub-steps
                                
                                //compute det(Rho(:,:,r)) [TODO: avoid re-computations]
-                               for (mwSize j=0; j<m; j++)
+                               for (int j=0; j<m; j++)
                                {
-                                       for (mwSize jj=0; jj<m; jj++)
-                                               matrixE->data[j*m+jj] = Rho[j*m*k+jj*k+r];
+                                       for (int jj=0; jj<m; jj++)
+                                               matrixE->data[j*m+jj] = Rho[ai(j,jj,r,m,m,k)];
                                }
                                gsl_linalg_LU_decomp(matrixE, permutation, &signum);
-                               Real detRhoR = gsl_linalg_LU_det(matrixE, signum);
-                               
+                               double detRhoR = gsl_linalg_LU_det(matrixE, signum);
+
                                //compute Y(i,:)*Rho(:,:,r)
-                               for (mwSize j=0; j<m; j++)
+                               for (int j=0; j<m; j++)
                                {
                                        YiRhoR[j] = 0.0;
-                                       for (mwSize u=0; u<m; u++)
-                                               YiRhoR[j] += Y[i*m+u] * Rho[u*m*k+j*k+r];
+                                       for (int u=0; u<m; u++)
+                                               YiRhoR[j] += Y[mi(i,u,n,m)] * Rho[ai(u,j,r,m,m,k)];
                                }
-                               
+
                                //compute X(i,:)*phi(:,:,r)
-                               for (mwSize j=0; j<m; j++)
+                               for (int j=0; j<m; j++)
                                {
                                        XiPhiR[j] = 0.0;
-                                       for (mwSize u=0; u<p; u++)
-                                               XiPhiR[j] += X[i*p+u] * phi[u*m*k+j*k+r];
+                                       for (int u=0; u<p; u++)
+                                               XiPhiR[j] += X[mi(i,u,n,p)] * phi[ai(u,j,r,p,m,k)];
                                }
-                               
+
                                //compute dotProduct < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) >
-                               Real dotProduct = 0.0;
-                               for (mwSize u=0; u<m; u++)
+                               double dotProduct = 0.0;
+                               for (int u=0; u<m; u++)
                                        dotProduct += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
-                               Real logGamIR = log(Pi[r]) + log(detRhoR) - 0.5*dotProduct;
-                               
+                               double logGamIR = log(Pi[r]) + log(detRhoR) - 0.5*dotProduct;
+
                                //Z(i) = index of max (gam(i,:))
                                if (logGamIR > maxLogGamIR)
                                {
@@ -240,23 +239,23 @@ void EMGrank(
                                }
                                sumLLF1 += exp(logGamIR) / pow(2*M_PI,m/2.0);
                        }
-                       
+
                        sumLogLLF2 += log(sumLLF1);
                }
-               
+
                // Assign output variable LLF
                *LLF = -invN * sumLogLLF2;
 
                //newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi))));
-               Real newDeltaPhi = 0.0;
-               for (mwSize j=0; j<p; j++)
+               double newDeltaPhi = 0.0;
+               for (int j=0; j<p; j++)
                {
-                       for (mwSize jj=0; jj<m; jj++)
+                       for (int jj=0; jj<m; jj++)
                        {
-                               for (mwSize r=0; r<k; r++)
+                               for (int r=0; r<k; r++)
                                {
-                                       Real tmpDist = fabs(phi[j*m*k+jj*k+r]-Phi[j*m*k+jj*k+r]) 
-                                               / (1.0+fabs(phi[j*m*k+jj*k+r]));
+                                       double tmpDist = fabs(phi[ai(j,jj,r,p,m,k)]-Phi[ai(j,jj,r,p,m,k)])
+                                               / (1.0+fabs(phi[ai(j,jj,r,p,m,k)]));
                                        if (tmpDist > newDeltaPhi)
                                                newDeltaPhi = tmpDist;
                                }
@@ -277,17 +276,17 @@ void EMGrank(
                sumDeltaPhi += newDeltaPhi;
 
                // update other local variables
-               for (mwSize j=0; j<m; j++)
+               for (int j=0; j<m; j++)
                {
-                       for (mwSize jj=0; jj<p; jj++)
+                       for (int jj=0; jj<p; jj++)
                        {
-                               for (mwSize r=0; r<k; r++)
-                                       Phi[j*m*k+jj*k+r] = phi[j*m*k+jj*k+r];
+                               for (int r=0; r<k; r++)
+                                       Phi[ai(j,jj,r,p,m,k)] = phi[ai(j,jj,r,p,m,k)];
                        }
                }
-        ite++;
+               ite++;
        }
-       
+
        //free memory
        free(hatBetaR);
        free(deltaPhi);
@@ -304,5 +303,5 @@ void EMGrank(
        free(Yr);
        free(tXrXr);
        free(tXrYr);
-    free(Z);
+       free(Z);
 }