-SelMix = setRefClass(
- Class = "selmix",
+Valse = setRefClass(
+ Class = "Valse",
fields = c(
# User defined
# regression data (size n*p, where n is the number of observations,
# and p is the number of regressors)
X = "numeric",
- # response data (size n*m, where n is the number of observations,
+ # response data (size n*m, where n is the number of observations,
# and m is the number of responses)
Y = "numeric",
#######################
initialize = function(X,Y,...)
{
- "Initialize SelMix object"
+ "Initialize Valse object"
callSuper(...)
- X <<- X;
- Y <<- Y;
+ X <<- X
+ Y <<- Y
gamma <<- ifelse (hasArg("gamma"), gamma, 1.)
mini <<- ifelse (hasArg("mini"), mini, as.integer(5))
maxi <<- ifelse (hasArg("maxi"), maxi, as.integer(10))
#smallEM initializes parameters by k-means and regression model in each component,
#doing this 20 times, and keeping the values maximizing the likelihood after 10
#iterations of the EM algorithm.
- init = initSmallEM(k,sx.X,sx.Y,sx.eps);
- phiInit <<- init$phi0;
- rhoInit <<- init$rho0;
- piInit <<- init$pi0;
- tauInit <<- init$tau0;
+ init = initSmallEM(k,X,Y,eps)
+ phiInit <<- init$phi0
+ rhoInit <<- init$rho0
+ piInit <<- init$pi0
+ tauInit <<- init$tau0
},
computeGridLambda = function()
"computation of the regularization grid"
#(according to explicit formula given by EM algorithm)
- gridLambda <<- gridLambda(sx.phiInit,sx.rhoInit,sx.piInit,sx.tauInit,sx.X,sx.Y,
- sx.gamma,sx.mini,sx.maxi,sx.eps);
+ gridLambda <<- gridLambda(phiInit,rhoInit,piInit,tauInit,X,Y,gamma,mini,maxi,eps)
},
computeRelevantParameters = function()
"Compute relevant parameters"
#select variables according to each regularization parameter
- #from the grid: sx.A1 corresponding to selected variables, and
- #sx.A2 corresponding to unselected variables.
- params = selectiontotale(sx.phiInit,sx.rhoInit,sx.piInit,sx.tauInit,sx.mini,sx.maxi,
- sx.gamma,sx.gridLambda,sx.X,sx.Y,sx.seuil,sx.eps);
+ #from the grid: A1 corresponding to selected variables, and
+ #A2 corresponding to unselected variables.
+ params = selectiontotale(
+ phiInit,rhoInit,piInit,tauInit,mini,maxi,gamma,gridLambda,X,Y,seuil,eps)
A1 <<- params$A1
A2 <<- params$A2
Rho <<- params$Rho
#compute parameter estimations, with the Maximum Likelihood
#Estimator, restricted on selected variables.
- res = constructionModelesLassoMLE(sx.phiInit,sx.rhoInit,sx.piInit,sx.tauInit,
- sx.mini,sx.maxi,sx.gamma,sx.gridLambda,sx.X,sx.Y,sx.seuil,sx.eps,sx.A1,sx.A2);
- return (list( phi=res$phi, rho=res$rho, pi=res$pi))
+ return ( constructionModelesLassoMLE(
+ phiInit,rhoInit,piInit,tauInit,mini,maxi,gamma,gridLambda,X,Y,seuil,eps,A1,A2) )
},
runProcedure2 = function()
#compute parameter estimations, with the Low Rank
#Estimator, restricted on selected variables.
- return (constructionModelesLassoRank(sx.Pi,sx.Rho,sx.mini,sx.maxi,sx.X,sx.Y,sx.eps,
- sx.A1,sx.rangmin,sx.rangmax)$phi)
+ return ( constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y,eps,
+ A1,rangmin,rangmax) )
},
- run = function(procedure)
+ run = function()
{
"main loop: over all k and all lambda"
# Run the all procedure, 1 with the
#maximum likelihood refitting, and 2 with the Low Rank refitting.
- p = dim(phiInit)[1]
- m = dim(phiInit)[2]
- for (k in kmin:kmax)
+ p = dim(phiInit)[1]
+ m = dim(phiInit)[2]
+ for (k in kmin:kmax)
+ {
+ print(k)
+ initParameters(k)
+ computeGridLambda()
+ computeRelevantParameters()
+ if (procedure == 1)
{
- print(k)
- initParameters(k)
- computeGridLambda()
- computeRelevantParameters()
- if (procedure == 1)
+ r1 = runProcedure1()
+ Phi2 = Phi
+ Rho2 = Rho
+ Pi2 = Pi
+ p = ncol(X)
+ m = ncol(Y)
+ if size(Phi2) == 0
+ {
+ Phi[,,1:k] = r1$phi
+ Rho[,,1:k] = r1$rho
+ Pi[1:k,] = r1$pi
+ } else
+ {
+ Phi = array(0., dim=c(p,m,kmax,dim(Phi2)[4]+dim(r1$phi)[4]))
+ Phi[,,1:(dim(Phi2)[3]),1:(dim(Phi2)[4])] = Phi2
+ Phi[,,1:k,dim(Phi2)[4]+1] = r1$phi
+ Rho = array(0., dim=c(m,m,kmax,dim(Rho2)[4]+dim(r1$rho)[4]))
+ Rho[,,1:(dim(Rho2)[3]),1:(dim(Rho2)[4])] = Rho2
+ Rho[,,1:k,dim(Rho2)[4]+1] = r1$rho
+ Pi = array(0., dim=c(kmax,dim(Pi2)[2]+dim(r1$pi)[2]))
+ Pi[1:nrow(Pi2),1:ncol(Pi2)] = Pi2
+ Pi[1:k,ncol(Pi2)+1] = r1$pi
+ }
+ } else
+ {
+ phi = runProcedure2()$phi
+ Phi2 = Phi
+ if (dim(Phi2)[1] == 0)
+ {
+ Phi(:,:,1:k,:) = phi
+ } else
{
- r1 = runProcedure1(sx)
- Phi2 = Phi
- Rho2 = Rho
- Pi2 = Pi
- p = ncol(X)
- m = ncol(Y)
- if size(Phi2) == 0 #TODO: continue translation MATLAB --> R
- Phi(:,:,1:k,:) = r1$phi;
- Rho(:,:,1:k,:) = r1$rho;
- Pi(1:k,:) = r1$pi;
- else
- Phi = zeros(p,m,sx.kmax,size(Phi2,4)+size(r1$phi,4));
- Phi(:,:,1:size(Phi2,3),1:size(Phi2,4)) = Phi2;
- Phi(:,:,1:k,size(Phi2,4)+1:end) = r1$phi;
- Rho = zeros(m,m,sx.kmax,size(Rho2,4)+size(r1$rho,4));
- Rho(:,:,1:size(Rho2,3),1:size(Rho2,4)) = Rho2;
- Rho(:,:,1:k,size(Rho2,4)+1:end) = r1$rho;
- Pi = zeros(sx.kmax,size(Pi2,2)+size(r1$pi,2));
- Pi(1:size(Pi2,1),1:size(Pi2,2)) = Pi2;
- Pi(1:k,size(Pi2,2)+1:end) = r1$pi;
- end
- else
- [phi] = runProcedure2(sx);
- phi
- Phi2 = sx.Phi;
- if size(Phi2,1) == 0
- sx.Phi(:,:,1:k,:) = phi;
- else
size(Phi2)
- sx.Phi = zeros(p,m,sx.kmax,size(Phi2,4)+size(phi,4));
- size(sx.Phi)
- sx.Phi(:,:,1:size(Phi2,3),1:size(Phi2,4)) = Phi2;
- sx.Phi(:,:,1:k,size(Phi2,4)+1:end) = phi;
- end
-
- end
-
-
- end
- end
-
+ Phi = zeros(p,m,kmax,size(Phi2,4)+size(phi,4))
+ size(Phi)
+ Phi(:,:,1:size(Phi2,3),1:size(Phi2,4)) = Phi2
+ Phi(:,:,1:k,size(Phi2,4)+1:end) = phi
+ }
+ }
+ }
+ }
+
##################################################
- #pruning: select only one (or a few best ?!) model
+ #TODO: pruning: select only one (or a few best ?!) model
##################################################
#
- # function[model] selectModel(sx)
+ # function[model] selectModel(
# #TODO
- # #model = sxModel(...);
+ # #model = odel(...)
# end
-
+
)
)
-PKG_CFLAGS=-g -I.
+#Debug flags
+PKG_CFLAGS=-g -I./sources
+
+#Prod flags:
+#PKG_CFLAGS=-O2 -I./sources
PKG_LIBS=-lm
#include <R.h>
#include <Rdefines.h>
-#include "sources/EMGLLF.h"
+#include "EMGLLF.h"
SEXP EMGLLF(
SEXP phiInit_,
) {
// Get matrices dimensions
int n = INTEGER(getAttrib(X_, R_DimSymbol))[0];
- SEXP dim = getAttrib(phiInit_, R_DimSymbol)
+ SEXP dim = getAttrib(phiInit_, R_DimSymbol);
int p = INTEGER(dim)[0];
int m = INTEGER(dim)[1];
int k = INTEGER(dim)[2];
PROTECT(pi = allocVector(REALSXP, k));
PROTECT(LLF = allocVector(REALSXP, maxi-mini+1));
PROTECT(S = allocArray(REALSXP, dimPhiS));
- double* pPhi=REAL(phi), pRho=REAL(rho), pPi=REAL(pi), pLLF=REAL(LLF), pS=REAL(S);
+ double *pPhi=REAL(phi), *pRho=REAL(rho), *pPi=REAL(pi), *pLLF=REAL(LLF), *pS=REAL(S);
////////////////////
// Call to EMGLLF //
#include <R.h>
#include <Rdefines.h>
-#include "sources/EMGLLF.h"
+#include "EMGrank.h"
SEXP EMGLLF(
SEXP Pi_,
SEXP dimX = getAttrib(X_, R_DimSymbol);
int n = INTEGER(dimX)[0];
int p = INTEGER(dimX)[1];
- SEXP dimRho = getAttrib(Rho_, R_DimSymbol)
+ SEXP dimRho = getAttrib(Rho_, R_DimSymbol);
int m = INTEGER(dimRho)[0];
int k = INTEGER(dimRho)[2];
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);
+ double *pPhi=REAL(phi), *pLLF=REAL(LLF);
/////////////////////
// Call to EMGrank //
#include <R.h>
#include <Rdefines.h>
-#include "sources/EMGLLF.h"
+#include "constructionModelesLassoMLE.h"
SEXP EMGLLF(
SEXP phiInit_,
) {
// Get matrices dimensions
int n = INTEGER(getAttrib(X_, R_DimSymbol))[0];
- SEXP dim = getAttrib(phiInit_, R_DimSymbol)
+ SEXP dim = getAttrib(phiInit_, R_DimSymbol);
int p = INTEGER(dim)[0];
int m = INTEGER(dim)[1];
int k = INTEGER(dim)[2];
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);
+ double *pPhi=REAL(phi), *pRho=REAL(rho), *pPi=REAL(pi), *pLvraisemblance=REAL(lvraisemblance);
/////////////////////////////////////////
// Call to constructionModelesLassoMLE //
#include <R.h>
#include <Rdefines.h>
-#include "sources/EMGLLF.h"
+#include "constructionModelesLassoRank.h"
SEXP EMGLLF(
SEXP Pi_,
SEXP tau_,
SEXP A1_,
SEXP rangmin_,
- SEXP rangmax
+ 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)
+ SEXP dimRho = getAttrib(Rho_, R_DimSymbol);
int m = INTEGER(dimRho)[0];
int k = INTEGER(dimRho)[2];
int L = INTEGER(getAttrib(A1_, R_DimSymbol))[1];
double* Rho = REAL(Rho_);
double* X = REAL(X_);
double* Y = REAL(Y_);
- double* A1 = REAL(A1_);
+ int* A1 = INTEGER(A1_);
/////////////
// OUTPUTS //
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);
+ double *pPhi=REAL(phi), *pLvraisemblance=REAL(lvraisemblance);
//////////////////////////////////////////
// Call to constructionModelesLassoRank //
#include <R.h>
#include <Rdefines.h>
-#include "sources/EMGLLF.h"
+#include "selectiontotale.h"
SEXP EMGLLF(
SEXP phiInit_,
SEXP dimX = getAttrib(X_, R_DimSymbol);
int n = INTEGER(dimX)[0];
int p = INTEGER(dimX)[1];
- SEXP dimRho = getAttrib(rhoInit_, R_DimSymbol)
+ SEXP dimRho = getAttrib(rhoInit_, R_DimSymbol);
int m = INTEGER(dimRho)[0];
int k = INTEGER(dimRho)[2];
int L = INTEGER(getAttrib(glambda_, R_LengthSymbol))[0];
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);
+ double *pA1=REAL(A1), *pA2=REAL(A2), *pRho=REAL(rho), *pPi=REAL(pi);
/////////////////////////////
// Call to selectiontotale //
-#include "EMGLLF.h"
+#include "utils.h"
+#include <stdlib.h>
#include <gsl/gsl_linalg.h>
// TODO: don't recompute indexes every time......
void EMGLLF(
// IN parameters
- const double* phiInit, // parametre initial de moyenne renormalisé
- const double* rhoInit, // parametre initial de variance renormalisé
- const double* piInit, // parametre initial des proportions
- const double* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon
- int mini, // nombre minimal d'itérations dans l'algorithme EM
- int maxi, // nombre maximal d'itérations dans l'algorithme EM
- double gamma, // valeur de gamma : puissance des proportions dans la pénalisation pour un Lasso adaptatif
+ const double* phiInit, // parametre initial de moyenne renormalisé
+ const double* rhoInit, // parametre initial de variance renormalisé
+ const double* piInit, // parametre initial des proportions
+ const double* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon
+ int mini, // nombre minimal d'itérations dans l'algorithme EM
+ int maxi, // nombre maximal d'itérations dans l'algorithme EM
+ double gamma, // puissance des proportions dans la pénalisation pour un Lasso adaptatif
double lambda, // valeur du paramètre de régularisation du Lasso
- const double* X, // régresseurs
- const double* Y, // réponse
- double tau, // seuil pour accepter la convergence
+ const double* X, // régresseurs
+ const double* Y, // réponse
+ double tau, // seuil pour accepter la convergence
// OUT parameters (all pointers, to be modified)
- double* phi, // parametre de moyenne renormalisé, calculé par l'EM
- double* rho, // parametre de variance renormalisé, calculé par l'EM
- double* pi, // parametre des proportions renormalisé, calculé par l'EM
- double* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres
- double* S,
+ double* phi, // parametre de moyenne renormalisé, calculé par l'EM
+ double* rho, // parametre de variance renormalisé, calculé par l'EM
+ double* pi, // parametre des proportions renormalisé, calculé par l'EM
+ double* LLF, // log vraisemblance associée à cet échantillon, pour les valeurs estimées des paramètres
+ double* S,
// additional size parameters
- int n, // nombre d'echantillons
- int p, // nombre de covariables
- int m, // taille de Y (multivarié)
- int k) // nombre de composantes dans le mélange
+ int n, // nombre d'echantillons
+ int p, // nombre de covariables
+ int m, // taille de Y (multivarié)
+ int k) // nombre de composantes dans le mélange
{
//Initialize outputs
copyArray(phiInit, phi, p*m*k);
double prodGam2logPi2 = 0.;
for (int v=0; v<k; v++)
prodGam2logPi2 += gam2[v] * log(pi2[v]);
- while (-invN*a + lambda*piPowGammaDotB < -invN*prodGam2logPi2 + lambda*pi2PowGammaDotB && kk<1000)
+ while (-invN*a + lambda*piPowGammaDotB < -invN*prodGam2logPi2 + lambda*pi2PowGammaDotB
+ && kk<1000)
{
//pi2=pi+0.1^kk*(1/n*gam2-pi);
for (int v=0; v<k; v++)
sumNy21 += nY21[ai(u,mm,r,n,m,k)];
nY2[mi(mm,r,m,k)] = sumNy21;
//rho(mm,mm,r)=((ps(mm,r)+sqrt(ps(mm,r)^2+4*nY2(mm,r)*(gam2(r))))/(2*nY2(mm,r)));
- rho[ai(mm,mm,k,m,m,k)] = ( ps[mi(mm,r,m,k)] + sqrt( ps[mi(mm,r,m,k)]*ps[mi(mm,r,m,k)]
+ rho[ai(mm,mm,k,m,m,k)] = ( ps[mi(mm,r,m,k)] + sqrt( ps[mi(mm,r,m,k)]*ps[mi(mm,r,m,k)]
+ 4*nY2[mi(mm,r,m,k)] * (gam2[r]) ) ) / (2*nY2[mi(mm,r,m,k)]);
}
}
{
for (int mm=0; mm<m; mm++)
{
- //sum(phi(1:j-1,mm,r).*transpose(Gram2(j,1:j-1,r)))+sum(phi(j+1:p,mm,r).*transpose(Gram2(j,j+1:p,r)))
+ //sum(phi(1:j-1,mm,r).*transpose(Gram2(j,1:j-1,r)))+sum(phi(j+1:p,mm,r)
+ // .*transpose(Gram2(j,j+1:p,r)))
double dotPhiGram2 = 0.0;
for (int u=0; u<j; u++)
dotPhiGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)];
for (int u=j+1; u<p; u++)
dotPhiGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)];
//S(j,r,mm)=-rho(mm,mm,r)*ps2(j,mm,r)+sum(phi(1:j-1,mm,r).*transpose(Gram2(j,1:j-1,r)))
- // +sum(phi(j+1:p,mm,r).*transpose(Gram2(j,j+1:p,r)));
+ // +sum(phi(j+1:p,mm,r).*transpose(Gram2(j,j+1:p,r)));
S[ai(j,mm,r,p,m,k)] = -rho[ai(mm,mm,r,m,m,k)] * ps2[ai(j,mm,r,p,m,k)] + dotPhiGram2;
if (fabs(S[ai(j,mm,r,p,m,k)]) <= n*lambda*pow(pi[r],gamma))
phi[ai(j,mm,r,p,m,k)] = 0;
else if (S[ai(j,mm,r,p,m,k)] > n*lambda*pow(pi[r],gamma))
- phi[ai(j,mm,r,p,m,k)] = (n*lambda*pow(pi[r],gamma) - S[ai(j,mm,r,p,m,k)])
+ phi[ai(j,mm,r,p,m,k)] = (n*lambda*pow(pi[r],gamma) - S[ai(j,mm,r,p,m,k)])
/ Gram2[ai(j,j,r,p,p,k)];
else
- phi[ai(j,mm,r,p,m,k)] = -(n*lambda*pow(pi[r],gamma) + S[ai(j,mm,r,p,m,k)])
+ phi[ai(j,mm,r,p,m,k)] = -(n*lambda*pow(pi[r],gamma) + S[ai(j,mm,r,p,m,k)])
/ Gram2[ai(j,j,r,p,p,k)];
}
}
{
//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 Y(i,:)*rho(:,:,r)
XiPhiR[u] += X[mi(i,v,n,p)] * phi[ai(v,u,r,p,m,k)];
}
- // compute dotProduct < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) >
+ //compute dotProduct
+ // < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) >
dotProducts[r] = 0.0;
for (int u=0; u<m; u++)
dotProducts[r] += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
-#ifndef select_EMGLLF_H
-#define select_EMGLLF_H
-
-#include "ioutils.h"
+#ifndef valse_EMGLLF_H
+#define valse_EMGLLF_H
void EMGLLF(
// IN parameters
- const Real* phiInit,
- const Real* rhoInit,
- const Real* piInit,
- const Real* gamInit,
- Int mini,
- Int maxi,
- Real gamma,
- Real lambda,
- const Real* X,
- const Real* Y,
- Real tau,
+ const double* phiInit,
+ const double* rhoInit,
+ const double* piInit,
+ const double* gamInit,
+ int mini,
+ int maxi,
+ double gamma,
+ double lambda,
+ const double* X,
+ const double* Y,
+ double tau,
// OUT parameters
- Real* phi,
- Real* rho,
- Real* pi,
- Real* LLF,
- Real* S,
+ double* phi,
+ double* rho,
+ double* pi,
+ double* LLF,
+ double* S,
// additional size parameters
- mwSize n,
- mwSize p,
- mwSize m,
- mwSize k);
+ int n,
+ int p,
+ int m,
+ int k);
#endif
-#include "EMGrank.h"
+#include <stdlib.h>
#include <gsl/gsl_linalg.h>
+#include "utils.h"
// Compute pseudo-inverse of a square matrix
static double* pinv(const double* matrix, int dim)
int* Z = (int*)calloc(n, sizeof(int));
//Initialize phi to zero, because some M loops might exit before phi affectation
- for (int i=0; i<p*m*k; i++)
- phi[i] = 0.0;
+ zeroArray(phi, p*m*k);
while (ite<mini || (ite<maxi && sumDeltaPhi>tau))
{
-#ifndef select_EMGrank_H
-#define select_EMGrank_H
-
-#include "ioutils.h"
+#ifndef valse_EMGrank_H
+#define valse_EMGrank_H
void EMGrank(
// IN parameters
- const Real* Pi,
- const Real* Rho,
- Int mini,
- Int maxi,
- const Real* X,
- const Real* Y,
- Real tau,
- const Int* rank,
+ const double* Pi,
+ const double* Rho,
+ int mini,
+ int maxi,
+ const double* X,
+ const double* Y,
+ double tau,
+ const int* rank,
// OUT parameters
- Real* phi,
- Real* LLF,
+ double* phi,
+ double* LLF,
// additional size parameters
- mwSize n,
- mwSize p,
- mwSize m,
- mwSize k);
+ int n,
+ int p,
+ int m,
+ int k);
#endif
#include "EMGLLF.h"
-#include "constructionModelesLassoMLE.h"
+#include "utils.h"
+#include <stdlib.h>
#include <gsl/gsl_linalg.h>
#include <omp.h>
-#include "omp_num_threads.h"
// TODO: comment on constructionModelesLassoMLE purpose
void constructionModelesLassoMLE(
-#ifndef select_constructionModelesLassoMLE_H
-#define select_constructionModelesLassoMLE_H
-
-#include "ioutils.h"
+#ifndef valse_constructionModelesLassoMLE_H
+#define valse_constructionModelesLassoMLE_H
void constructionModelesLassoMLE(
- // IN parameters
- const Real* phiInit,
- const Real* rhoInit,
- const Real* piInit,
- const Real* gamInit,
- Int mini,
- Int maxi,
- Real gamma,
- const Real* glambda,
- const Real* X,
- const Real* Y,
- Real seuil,
- Real tau,
- const Int* A1,
- const Int* A2,
+ // IN parameters
+ const double* phiInit,
+ const double* rhoInit,
+ const double* piInit,
+ const double* gamInit,
+ int mini,
+ int maxi,
+ double gamma,
+ const double* glambda,
+ const double* X,
+ const double* Y,
+ double seuil,
+ double tau,
+ const int* A1,
+ const int* A2,
// OUT parameters
- Real* phi,
- Real* rho,
- Real* pi,
- Real* lvraisemblance,
+ double* phi,
+ double* rho,
+ double* pi,
+ double* lvraisemblance,
// additional size parameters
- mwSize n,
- mwSize p,
- mwSize m,
- mwSize k,
- mwSize L);
+ int n,
+ int p,
+ int m,
+ int k,
+ int L);
#endif
-#include "EMGrank.h"
-#include "constructionModelesLassoRank.h"
-#include <gsl/gsl_linalg.h>
+#include <stdlib.h>
#include <omp.h>
-#include "omp_num_threads.h"
+#include <gsl/gsl_linalg.h>
+#include "EMGrank.h"
+#include "utils.h"
// TODO: comment on constructionModelesLassoRank purpose
void constructionModelesLassoRank(
-#ifndef select_constructionModelesLassoRank_H
-#define select_constructionModelesLassoRank_H
-
-#include "ioutils.h"
+#ifndef valse_constructionModelesLassoRank_H
+#define valse_constructionModelesLassoRank_H
// Main job on raw inputs (after transformation from mxArray)
void constructionModelesLassoRank(
- // IN parameters
- const Real* Pi,
- const Real* Rho,
- Int mini,
- Int maxi,
- const Real* X,
- const Real* Y,
- Real tau,
- const Int* A1,
- Int rangmin,
- Int rangmax,
+ // IN parameters
+ const double* Pi,
+ const double* Rho,
+ int mini,
+ int maxi,
+ const double* X,
+ const double* Y,
+ double tau,
+ const int* A1,
+ int rangmin,
+ int rangmax,
// OUT parameters
- Real* phi,
- Real* lvraisemblance,
- // additional size parameters
- mwSize n,
- mwSize p,
- mwSize m,
- mwSize k,
- mwSize L);
+ double* phi,
+ double* lvraisemblance,
+ // additional size parameters
+ int n,
+ int p,
+ int m,
+ int k,
+ int L);
#endif
-#include "selectiontotale.h"
-#include "EMGLLF.h"
+#include <stdlib.h>
#include <omp.h>
-#include "omp_num_threads.h"
+#include "EMGLLF.h"
+#include "utils.h"
// Main job on raw inputs (after transformation from mxArray)
void selectiontotale(
- // IN parameters
- const Real* phiInit, // parametre initial de moyenne renormalisé
- const Real* rhoInit, // parametre initial de variance renormalisé
- const Real* piInit, // parametre initial des proportions
- const Real* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon
- Int mini, // nombre minimal d'itérations dans lambdaIndex'algorithme EM
- Int maxi, // nombre maximal d'itérations dans lambdaIndex'algorithme EM
- Real gamma, // valeur de gamma : puissance des proportions dans la pénalisation pour un Lasso adaptatif
- const Real* glambda, // valeur des paramètres de régularisation du Lasso
- const Real* X, // régresseurs
- const Real* Y, // réponse
- Real seuil, // seuil pour prendre en compte une variable
- Real tau, // seuil pour accepter la convergence
+ // IN parameters
+ const double* phiInit, // parametre initial de moyenne renormalisé
+ const double* rhoInit, // parametre initial de variance renormalisé
+ const double* piInit,// parametre initial des proportions
+ const double* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon
+ int mini, // nombre minimal d'itérations dans lambdaIndex'algorithme EM
+ int maxi, // nombre maximal d'itérations dans lambdaIndex'algorithme EM
+ double gamma, // valeur de gamma : puissance des proportions dans la pénalisation pour un Lasso adaptatif
+ const double* glambda, // valeur des paramètres de régularisation du Lasso
+ const double* X,// régresseurs
+ const double* Y,// réponse
+ double seuil, // seuil pour prendre en compte une variable
+ double tau, // seuil pour accepter la convergence
// OUT parameters (all pointers, to be modified)
- Int* A1, // matrice des coefficients des parametres selectionnes
- Int* A2, // matrice des coefficients des parametres non selectionnes
- Real* Rho, // estimateur ainsi calculé par le Lasso
- Real* Pi, // estimateur ainsi calculé par le Lasso
+ int* A1, // matrice des coefficients des parametres selectionnes
+ int* A2, // matrice des coefficients des parametres non selectionnes
+ double* Rho,// estimateur ainsi calculé par le Lasso
+ double* Pi,// estimateur ainsi calculé par le Lasso
// additional size parameters
- mwSize n, // taille de lambdaIndex'echantillon
- mwSize p, // nombre de covariables
- mwSize m, // taille de Y (multivarié)
- mwSize k, // nombre de composantes
- mwSize L) // taille de glambda
+ int n,// taille de lambdaIndex'echantillon
+ int p,// nombre de covariables
+ int m,// taille de Y (multivarié)
+ int k,// nombre de composantes
+ int L) // taille de glambda
{
// Fill outputs with zeros: they might not be assigned
for (int u=0; u<p*(m+1)*L; u++)
Rho[u] = 0.0;
for (int u=0; u<k*L; u++)
Pi[u] = 0.0;
-
+
//initiate parallel section
- mwSize lambdaIndex;
+ int lambdaIndex;
omp_set_num_threads(OMP_NUM_THREADS);
#pragma omp parallel default(shared) private(lambdaIndex)
{
for (lambdaIndex=0; lambdaIndex<L; lambdaIndex++)
{
//allocate output variables
- Real* phi = (Real*)malloc(p*m*k*sizeof(Real));
- Real* rho = (Real*)malloc(m*m*k*sizeof(Real));
- Real* pi = (Real*)malloc(k*sizeof(Real));
- Real* LLF = (Real*)malloc(maxi*sizeof(Real));
- Real* S = (Real*)malloc(p*m*k*sizeof(Real));
+ double* phi = (double*)malloc(p*m*k*sizeof(double));
+ double* rho = (double*)malloc(m*m*k*sizeof(double));
+ double* pi = (double*)malloc(k*sizeof(double));
+ double* LLF = (double*)malloc(maxi*sizeof(double));
+ double* S = (double*)malloc(p*m*k*sizeof(double));
EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda[lambdaIndex],X,Y,tau,
phi,rho,pi,LLF,S,
n,p,m,k);
free(LLF);
free(S);
-
+
//Si un des coefficients est supérieur au seuil, on garde cette variable
- mwSize* selectedVariables = (mwSize*)calloc(p*m,sizeof(mwSize));
- mwSize* discardedVariables = (mwSize*)calloc(p*m,sizeof(mwSize));
+ int* selectedVariables = (int*)calloc(p*m,sizeof(int));
+ int* discardedVariables = (int*)calloc(p*m,sizeof(int));
int atLeastOneSelectedVariable = 0;
- for (
mwSize j=0; j<p; j++)
+ for (
int j=0; j<p; j++)
{
- mwSize cpt = 0;
- mwSize cpt2 = 0;
- for (mwSize jj=0; jj<m; jj++)
+ int cpt = 0;
+ int cpt2 = 0;
+ for (int jj=0; jj<m; jj++)
{
- Real maxPhi = 0.0;
- for (mwSize r=0; r<k; r++)
+ double maxPhi = 0.0;
+ for (int r=0; r<k; r++)
{
- if (fabs(phi[j*m*k+jj*k+r]) > maxPhi)
- maxPhi = fabs(phi[j*m*k+jj*k+r]);
+ if (fabs(phi[ai(j,jj,r,p,m,k)]) > maxPhi)
+ maxPhi = fabs(phi[ai(j,jj,r,p,m,k)]);
}
if (maxPhi > seuil)
{
- selectedVariables[j*m+cpt] = jj+1;
+ selectedVariables[mi(j,cpt,p,m)] = jj+1;
atLeastOneSelectedVariable = 1;
cpt++;
}
else
{
- discardedVariables[j*m+cpt2] = jj+1;
+ discardedVariables[mi(j,cpt2,p,m)] = jj+1;
cpt2++;
}
}
}
free(phi);
-
+
if (atLeastOneSelectedVariable)
{
- Int* vec = (Int*)malloc(p*sizeof(Int));
- mwSize vecSize = 0;
- for (
mwSize j=0; j<p; j++)
+ int* vec = (int*)malloc(p*sizeof(int));
+ int vecSize = 0;
+ for (
int j=0; j<p; j++)
{
- if (selectedVariables[j*m+0] != 0)
+ if (selectedVariables[mi(j,0,p,m)] != 0)
vec[vecSize++] = j;
}
//A1
- for (
mwSize j=0; j<p; j++)
- A1[j*(m+1)*L+0*L+lambdaIndex] = (j < vecSize ? vec[j]+1 : 0);
- for (mwSize j=0; j<vecSize; j++)
+ for (
int j=0; j<p; j++)
+ A1[ai(j,0,lambdaIndex,p,m+1,L)] = (j < vecSize ? vec[j]+1 : 0);
+ for (int j=0; j<vecSize; j++)
{
- for (mwSize jj=1; jj<=m; jj++)
- A1[j*(m+1)*L+jj*L+lambdaIndex] = selectedVariables[vec[j]*m+jj-1];
+ for (int jj=1; jj<=m; jj++)
+ A1[ai(j,jj,lambdaIndex,p,m+1,L)] = selectedVariables[mi(vec[j],jj-1,p,m)];
}
//A2
- for (
mwSize j=0; j<p; j++)
- A2[j*(m+1)*L+0*L+lambdaIndex] = j+1;
- for (
mwSize j=0; j<p; j++)
+ for (
int j=0; j<p; j++)
+ A2[ai(j,0,lambdaIndex,p,m+1,L)] = j+1;
+ for (
int j=0; j<p; j++)
{
- for (mwSize jj=1; jj<=m; jj++)
- A2[j*(m+1)*L+jj*L+lambdaIndex] = discardedVariables[j*m+jj-1];
+ for (int jj=1; jj<=m; jj++)
+ A2[ai(j,jj,lambdaIndex,p,m+1,L)] = discardedVariables[mi(j,jj-1,p,m)];
}
//Rho
- for (mwSize j=0; j<m; j++)
+ for (int j=0; j<m; j++)
{
- for (mwSize jj=0; jj<m; jj++)
+ for (int jj=0; jj<m; jj++)
{
- for (mwSize r=0; r<k; r++)
- Rho[j*m*k*L+jj*k*L+r*L+lambdaIndex] = rho[j*m*k+jj*k+r];
+ for (int r=0; r<k; r++)
+ Rho[ai4(j,jj,r,lambdaIndex,m,m,k,L)] = rho[ai(j,jj,r,m,m,k)];
}
}
//Pi
- for (mwSize r=0; r<k; r++)
- Pi[r*L+lambdaIndex] = pi[r];
+ for (int r=0; r<k; r++)
+ Pi[mi(r,lambdaIndex,k,L)] = pi[r];
free(vec);
}
-
+
free(selectedVariables);
free(discardedVariables);
free(rho);
-#ifndef select_selectiontotale_H
-#define select_selectiontotale_H
-
-#include "ioutils.h"
+#ifndef valse_selectiontotale_H
+#define valse_selectiontotale_H
// Main job on raw inputs (after transformation from mxArray)
void selectiontotale(
- // IN parameters
- const Real* phiInit,
- const Real* rhoInit,
- const Real* piInit,
- const Real* gamInit,
- Int mini,
- Int maxi,
- Real gamma,
- const Real* glambda,
- const Real* X,
- const Real* Y,
- Real seuil,
- Real tau,
+ // IN parameters
+ const double* phiInit,
+ const double* rhoInit,
+ const double* piInit,
+ const double* gamInit,
+ int mini,
+ int maxi,
+ double gamma,
+ const double* glambda,
+ const double* X,
+ const double* Y,
+ double seuil,
+ double tau,
// OUT parameters
- Int* A1,
- Int* A2,
- Real* Rho,
- Real* Pi,
+ int* A1,
+ int* A2,
+ double* Rho,
+ double* Pi,
// additional size parameters
- mwSize n,
- mwSize p,
- mwSize m,
- mwSize k,
- mwSize L);
-
+ int n,
+ int p,
+ int m,
+ int k,
+ int L);
+
#endif
--- /dev/null
+#ifndef valse_utils_H
+#define valse_utils_H
+
+/*******************
+ * tune parallelism
+ *******************/
+
+// Number of OpenMP threads
+#define OMP_NUM_THREADS 8
+
+// CHUNK_SIZE = number of lambda values to be treated sequentially by a single core
+#define CHUNK_SIZE 1
+
+/*******************************
+ * Matrix and arrays indexation
+ *******************************/
+
+// Matrix Index ; TODO? ncol unused
+#define mi(i,j,nrow,ncol)\
+ j*nrow + i
+
+// Array Index ; TODO? d3 unused
+#define ai(i,j,k,d1,d2,d3)\
+ k*d1*d2 + j*d1 + i
+
+// Array4 Index ; TODO? ...
+#define ai4(i,j,k,m,d1,d2,d3,d4)\
+ m*d1*d2*d3 + k*d1*d2 + j*d1 + i
+
+// Array5 Index ; TODO? ...
+#define ai5(i,j,k,m,n,d1,d2,d3,d4,d5)\
+ n*d1*d2*d3*d4 + m*d1*d2*d3 + k*d1*d2 + j*d1 + i
+
+/*************************
+ * Array copy & "zeroing"
+ ************************/
+
+// Fill an array with zeros
+#define zeroArray(array, size)\
+{\
+ for (int u=0; u<size; u++)\
+ array[u] = 0;\
+}
+
+// Copy an 1D array
+#define copyArray(array, copy, size)\
+{\
+ for (int u=0; u<size; u++)\
+ copy[u] = array[u];\
+}
+
+#endif
+++ /dev/null
-#ifndef select_ioutils_H
-#define select_ioutils_H
-
-// Fill an array with zeros
-#define zeroArray(array, size)\
-{\
- for (int u=0; u<size; u++)\
- array[u] = 0;\
-}
-
-// Copy an 1D array
-#define copyArray(array, copy, size)\
-{\
- for (int u=0; u<size; u++)\
- copy[u] = array[u];\
-}
-
-// Matrix Index ; TODO? ncol unused
-#define mi(i,j,nrow,ncol)\
- j*nrow + i
-
-// Array Index ; TODO? d3 unused
-#define ai(i,j,k,d1,d2,d3)\
- k*d1*d2 + j*d1 + i
-
-#endif
+++ /dev/null
-#ifndef tune_parallelism_H
-#define tune_parallelism_H
-
-// Number of OpenMP threads
-#define OMP_NUM_THREADS 8
-
-// CHUNK_SIZE = number of lambda values to be treated sequentially by a single core
-#define CHUNK_SIZE 1
-
-#endif
projects / valse.git / commitdiff