From 552b00e200e8a990c1247989d29e98d4ae8679f3 Mon Sep 17 00:00:00 2001 From: Benjamin Auder Date: Thu, 15 Dec 2016 13:17:04 +0100 Subject: [PATCH] finish adapters ; almost finished sources --- src/adapters/a.EMGrank.c | 120 ++++----- src/adapters/a.constructionModelesLassoMLE.c | 196 ++++++-------- src/adapters/a.constructionModelesLassoRank.c | 146 ++++------- src/adapters/a.selectiontotale.c | 178 +++++-------- src/sources/EMGrank.c | 247 +++++++++--------- 5 files changed, 374 insertions(+), 513 deletions(-) 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 - -// 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 +#include +#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 - -#include - -// 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 +#include +#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 +#include +#include +#include "sources/EMGLLF.h" -#include - -// 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 - -// 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 +#include +#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 // 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; idata[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; idata[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; itau)) - { + { ///////////// // Etape M // ///////////// //M step: Mise à jour de Beta (et donc phi) - for (mwSize r=0; rdata[j*m+jj] = dotProduct; + double dotProduct = 0.0; + for (int u=0; udata[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]; jdata[j] = 0.0; //[intermediate step] Compute hatBetaR = U * S * t(V) - Real* U = matrixM->data; - for (mwSize j=0; jdata; + for (int j=0; jdata[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; jdata[j*m+jj] = Rho[j*m*k+jj*k+r]; + for (int jj=0; jjdata[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 - Real dotProduct = 0.0; - for (mwSize u=0; u 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 newDeltaPhi) newDeltaPhi = tmpDist; } @@ -277,17 +276,17 @@ void EMGrank( sumDeltaPhi += newDeltaPhi; // update other local variables - for (mwSize j=0; j