-#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;
}
-#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;
}
-#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;
}
-#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;
}
#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;
}
// 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);
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)
{
}
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;
}
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);
free(Yr);
free(tXrXr);
free(tXrYr);
- free(Z);
+ free(Z);
}
projects / valse.git / commitdiff