From: Benjamin Auder <benjamin.auder@somewhere>
Date: Thu, 15 Dec 2016 12:17:04 +0000 (+0100)
Subject: finish adapters ; almost finished sources
X-Git-Url: https://git.auder.net/variants/Chakart/pieces/current/css/doc/R.css?a=commitdiff_plain;h=552b00e200e8a990c1247989d29e98d4ae8679f3;p=valse.git

finish adapters ; almost finished sources
---

diff --git a/src/adapters/a.EMGrank.c b/src/adapters/a.EMGrank.c
index 778ef67..f7fc192 100644
--- a/src/adapters/a.EMGrank.c
+++ b/src/adapters/a.EMGrank.c
@@ -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;
 }
diff --git a/src/adapters/a.constructionModelesLassoMLE.c b/src/adapters/a.constructionModelesLassoMLE.c
index a331c67..7578664 100644
--- a/src/adapters/a.constructionModelesLassoMLE.c
+++ b/src/adapters/a.constructionModelesLassoMLE.c
@@ -1,142 +1,92 @@
-#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;
 }
diff --git a/src/adapters/a.constructionModelesLassoRank.c b/src/adapters/a.constructionModelesLassoRank.c
index 96908c2..359290c 100644
--- a/src/adapters/a.constructionModelesLassoRank.c
+++ b/src/adapters/a.constructionModelesLassoRank.c
@@ -1,109 +1,79 @@
-#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;
 }
diff --git a/src/adapters/a.selectiontotale.c b/src/adapters/a.selectiontotale.c
index b4b8f4c..1124d71 100644
--- a/src/adapters/a.selectiontotale.c
+++ b/src/adapters/a.selectiontotale.c
@@ -1,127 +1,87 @@
-#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;
 }
diff --git a/src/sources/EMGrank.c b/src/sources/EMGrank.c
index 721c0c2..2e85ce6 100644
--- a/src/sources/EMGrank.c
+++ b/src/sources/EMGrank.c
@@ -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);
 }