# Optionally user defined (some default values)
# power in the penalty
- gamma = "double",
+ gamma = "numeric",
# minimum number of iterations for EM algorithm
mini = "integer",
# maximum number of iterations for EM algorithm
maxi = "integer",
# threshold for stopping EM algorithm
- eps = "double",
+ eps = "numeric",
# minimum number of components in the mixture
kmin = "integer",
# maximum number of components in the mixture
# Computed through the workflow
# initialisation for the reparametrized conditional mean parameter
- phiInit,
+ phiInit = "numeric",
# initialisation for the reparametrized variance parameter
- rhoInit,
+ rhoInit = "numeric",
# initialisation for the proportions
- piInit,
+ piInit = "numeric",
# initialisation for the allocations probabilities in each component
- tauInit,
+ tauInit = "numeric",
# values for the regularization parameter grid
- gridLambda = c(),
+ gridLambda = "numeric",
# je ne crois pas vraiment qu'il faille les mettre en sortie, d'autant plus qu'on construit
# une matrice A1 et A2 pour chaque k, et elles sont grandes, donc ca coute un peu cher ...
- A1,
- A2,
+ A1 = "integer",
+ A2 = "integer",
# collection of estimations for the reparametrized conditional mean parameters
- Phi,
+ Phi = "numeric",
# collection of estimations for the reparametrized variance parameters
- Rho,
+ Rho = "numeric",
# collection of estimations for the proportions parameters
- Pi,
+ Pi = "numeric",
- #immutable
- seuil = 1e-15
+ #immutable (TODO:?)
+ seuil = "numeric"
),
methods = list(
kmax <<- ifelse (hasArg("kmax"), kmax, as.integer(3))
rangmin <<- ifelse (hasArg("rangmin"), rangmin, as.integer(2))
rangmax <<- ifelse (hasArg("rangmax"), rangmax, as.integer(3))
+ seuil <<- 1e-15 #immutable (TODO:?)
},
##################################
Pi2 = Pi
p = ncol(X)
m = ncol(Y)
- if size(Phi2) == 0
+ if (is.null(dim(Phi2))) #test was: size(Phi2) == 0
{
- Phi[,,1:k] = r1$phi
- Rho[,,1:k] = r1$rho
- Pi[1:k,] = r1$pi
+ 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
+ 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
{
Phi2 = Phi
if (dim(Phi2)[1] == 0)
{
- Phi(:,:,1:k,:) = phi
+ Phi[,,1:k,] <<- phi
} else
{
- size(Phi2)
- 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
+ Phi <<- array(0., dim=c(p,m,kmax,dim(Phi2)[4]+dim(phi)[4]))
+ Phi[,,1:(dim(Phi2)[3]),1:(dim(Phi2)[4])] <<- Phi2
+ Phi[,,1:k,-(1:(dim(Phi2)[4]))] <<- phi
}
}
}
-#' selectVaribles
-#' It is a function which construct, for a given lambda, the sets of
+#' selectVaribles
+#' It is a function which construct, for a given lambda, the sets of
#' relevant variables and irrelevant variables.
#'
#' @param phiInit an initial estimator for phi (size: p*m*k)
#' @param rhoInit an initial estimator for rho (size: m*m*k)
-#' @param piInit an initial estimator for pi (size : k)
+#' @param piInit an initial estimator for pi (size : k)
#' @param gamInit an initial estimator for gamma
-#' @param mini minimum number of iterations in EM algorithm
-#' @param maxi maximum number of iterations in EM algorithm
-#' @param gamma power in the penalty
+#' @param mini minimum number of iterations in EM algorithm
+#' @param maxi maximum number of iterations in EM algorithm
+#' @param gamma power in the penalty
#' @param glambda grid of regularization parameters
-#' @param X matrix of regressors
-#' @param Y matrix of responses
-#' @param thres threshold to consider a coefficient to be equal to 0
-#' @param tau threshold to say that EM algorithm has converged
+#' @param X matrix of regressors
+#' @param Y matrix of responses
+#' @param thres threshold to consider a coefficient to be equal to 0
+#' @param tau threshold to say that EM algorithm has converged
#'
#' @return
#' @export
#'
#' @examples
selectVariables <- function(phiInit,rhoInit,piInit,gamInit,
- mini,maxi,gamma,glambda,X,Y,thres,tau){
-
- dimphi <- dim(phiInit)
- p <- dimPhi[1]
- m <- dimPhi[2]
- k <- dimPhi[3]
- L <- length(glambda);
- A1 <- array(0, dim <- c(p,m+1,L))
- A2 <- array(0, dim <- c(p,m+1,L))
- Rho <- array(0, dim <- c(m,m,k,L))
- Pi <- array(0, dim <- c(k,L));
-
- # For every lambda in gridLambda, comutation of the coefficients
- for (lambdaIndex in c(1:L)) {
- Res <- EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,
- gamma,glambda[lambdaIndex],X,Y,tau);
- phi <- Res$phi
- rho <- Res$rho
- pi <- Res$pi
-
- # If a coefficient is larger than the threshold, we keep it
- selectedVariables <- array(0, dim = c(p,m))
- discardedVariables <- array(0, dim = c(p,m))
- atLeastOneSelectedVariable <- false
- for (j in c(1:p)){
- cpt <- 1
- cpt2 <-1
- for (mm in c(1:m)){
- if (max(abs(phi[j,mm,])) > thres){
- selectedVariables[j,cpt] <- mm
- cpt <- cpt+1
- atLeastOneSelectedVariable <- true
- } else{
- discardedVariables[j,cpt2] <- mm
- cpt2 <- cpt2+1
- }
- }
- }
-
- # If no coefficients have been selected, we provide the zero matrix
- # We delete zero coefficients: vec = indices of zero values
- if atLeastOneSelectedVariable{
- vec <- c()
- for (j in c(1:p)){
- if (selectedVariables(j,1) =! 0){
- vec <- c(vec,j)
- }
- }
- # Else, we provide the indices of relevant coefficients
- A1[,1,lambdaIndex] <- c(vec,rep(0,p-length(vec)))
- A1[1:length(vec),2:(m+1),lambdaIndex] <- selectedVariables[vec,]
- A2[,1,lambdaIndex] <- 1:p
- A2[,2:(m+1),lambdaIndex] <- discardedVariables
- Rho[,,,lambdaIndex] <- rho
- Pi[,lambdaIndex] <- pi
- }
-
- }
- return(res = list(A1 = A1, A2 = A2 , Rho = Rho, Pi = Pi))
-}
\ No newline at end of file
+ mini,maxi,gamma,glambda,X,Y,thres,tau)
+{
+ dimphi <- dim(phiInit)
+ p <- dimPhi[1]
+ m <- dimPhi[2]
+ k <- dimPhi[3]
+ L <- length(glambda);
+ A1 <- array(0, dim <- c(p,m+1,L))
+ A2 <- array(0, dim <- c(p,m+1,L))
+ Rho <- array(0, dim <- c(m,m,k,L))
+ Pi <- array(0, dim <- c(k,L));
+
+ # For every lambda in gridLambda, comutation of the coefficients
+ for (lambdaIndex in c(1:L))
+ {
+ Res <- EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,
+ gamma,glambda[lambdaIndex],X,Y,tau);
+ phi <- Res$phi
+ rho <- Res$rho
+ pi <- Res$pi
+
+ # If a coefficient is larger than the threshold, we keep it
+ selectedVariables <- array(0, dim = c(p,m))
+ discardedVariables <- array(0, dim = c(p,m))
+ atLeastOneSelectedVariable <- false
+ for (j in c(1:p))
+ {
+ cpt <- 1
+ cpt2 <-1
+ for (mm in c(1:m))
+ {
+ if (max(abs(phi[j,mm,])) > thres)
+ {
+ selectedVariables[j,cpt] <- mm
+ cpt <- cpt+1
+ atLeastOneSelectedVariable <- true
+ } else
+ {
+ discardedVariables[j,cpt2] <- mm
+ cpt2 <- cpt2+1
+ }
+ }
+ }
+
+ # If no coefficients have been selected, we provide the zero matrix
+ # We delete zero coefficients: vec = indices of zero values
+ if (atLeastOneSelectedVariable)
+ {
+ vec <- c()
+ for (j in c(1:p))
+ {
+ if (selectedVariables(j,1) != 0)
+ vec <- c(vec,j)
+ # Else ( NOTE: [auder] else ?! TODO: explain? )
+ # we provide the indices of relevant coefficients
+ A1[,1,lambdaIndex] <- c(vec,rep(0,p-length(vec)))
+ A1[1:length(vec),2:(m+1),lambdaIndex] <- selectedVariables[vec,]
+ A2[,1,lambdaIndex] <- 1:p
+ A2[,2:(m+1),lambdaIndex] <- discardedVariables
+ Rho[,,,lambdaIndex] <- rho
+ Pi[,lambdaIndex] <- pi
+ }
+ }
+ }
+
+ return(res = list(A1 = A1, A2 = A2 , Rho = Rho, Pi = Pi))
+}
// Call to EMGLLF //
////////////////////
- EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau,
+ EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau,
pPhi,pRho,pPi,pLLF,pS,
n,p,m,k);
for (int i=0; i<5; 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, 0, phi);
+ SET_VECTOR_ELT(listParams, 1, rho);
+ SET_VECTOR_ELT(listParams, 2, pi);
SET_VECTOR_ELT(listParams, 3, LLF);
- SET_ARRAY_ELT(listParams, 4, S);
+ SET_VECTOR_ELT(listParams, 4, S);
UNPROTECT(9);
return listParams;
#include <Rdefines.h>
#include "EMGrank.h"
-SEXP EMGLLF(
+SEXP EMGrank(
SEXP Pi_,
SEXP Rho_,
SEXP mini_,
double* Rho = REAL(Rho_);
double* X = REAL(X_);
double* Y = REAL(Y_);
- double* rank = REAL(rank_);
+ int* rank = INTEGER(rank_);
/////////////
// OUTPUTS //
/////////////
- SEXP phi, LLF;
+ SEXP phi, LLF, dimPhi;
PROTECT(dimPhi = allocVector(INTSXP, 3));
int* pDimPhi = INTEGER(dimPhi);
pDimPhi[0] = p; pDimPhi[1] = m; pDimPhi[2] = k;
// Call to EMGrank //
/////////////////////
- EMGrank(Pi, Rho, mini, maxi, X, Y, tau, rank,
+ EMGrank_core(Pi, Rho, mini, maxi, X, Y, tau, rank,
pPhi,pLLF,
n,p,m,k);
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, 0, phi);
SET_VECTOR_ELT(listParams, 1, LLF);
UNPROTECT(5);
#include <Rdefines.h>
#include "constructionModelesLassoMLE.h"
-SEXP EMGLLF(
+SEXP constructionModelesLassoMLE(
SEXP phiInit_,
SEXP rhoInit_,
SEXP piInit_,
int p = INTEGER(dim)[0];
int m = INTEGER(dim)[1];
int k = INTEGER(dim)[2];
- int L = INTEGER(getAttrib(glambda_, R_LengthSymbol))[0];
+ int L = length(glambda_);
////////////
// INPUTS //
double* glambda = REAL(glambda_);
double* X = REAL(X_);
double* Y = REAL(Y_);
- double* A1 = REAL(A1_);
- double* A2 = REAL(A2_);
+ int* A1 = INTEGER(A1_);
+ int* A2 = INTEGER(A2_);
/////////////
// OUTPUTS //
// Call to constructionModelesLassoMLE //
/////////////////////////////////////////
- constructionModelesLassoMLE(
+ constructionModelesLassoMLE_core(
phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2,
pPhi,pRho,pPi,pLvraisemblance,
n,p,m,k,L);
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, 0, phi);
+ SET_VECTOR_ELT(listParams, 1, rho);
+ SET_VECTOR_ELT(listParams, 2, pi);
SET_VECTOR_ELT(listParams, 3, lvraisemblance);
UNPROTECT(8);
#include <Rdefines.h>
#include "constructionModelesLassoRank.h"
-SEXP EMGLLF(
+SEXP constructionModelesLassoRank(
SEXP Pi_,
SEXP Rho_,
SEXP mini_,
// Call to constructionModelesLassoRank //
//////////////////////////////////////////
- constructionModelesLassoRank(
+ constructionModelesLassoRank_core(
Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax,
pPhi,pLvraisemblance,
n,p,m,k,L);
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, 0, phi);
SET_VECTOR_ELT(listParams, 1, lvraisemblance);
UNPROTECT(5);
#include <Rdefines.h>
#include "selectiontotale.h"
-SEXP EMGLLF(
+SEXP selectiontotale(
SEXP phiInit_,
SEXP rhoInit_,
SEXP piInit_,
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];
+ SEXP dimRhoInit = getAttrib(rhoInit_, R_DimSymbol);
+ int m = INTEGER(dimRhoInit)[0];
+ int k = INTEGER(dimRhoInit)[2];
+ int L = length(glambda_);
////////////
// INPUTS //
double tau = NUMERIC_VALUE(tau_);
// Get pointers from SEXP arrays ; WARNING: by columns !
- double* piInit = REAL(phiInit_);
+ double* phiInit = REAL(phiInit_);
double* rhoInit = REAL(rhoInit_);
double* piInit = REAL(piInit_);
double* gamInit = REAL(gamInit_);
// OUTPUTS //
/////////////
- int Size = pow(rangmax-rangmin+1,k);
- SEXP A1, A2, rho, pi, dimA;
+ SEXP A1, A2, rho, pi, dimA, dimRho;
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(A1 = allocArray(INTSXP, dimA));
+ PROTECT(A2 = allocArray(INTSXP, dimA));
+ PROTECT(dimRho = allocVector(INTSXP, 4));
+ int* pDimRho = INTEGER(dimRho);
+ pDimRho[0] = m; pDimRho[1] = m; pDimRho[2] = k; pDimRho[3] = L;
+ PROTECT(rho = allocArray(REALSXP, dimRho));
PROTECT(pi = allocMatrix(REALSXP, k, L));
- double *pA1=REAL(A1), *pA2=REAL(A2), *pRho=REAL(rho), *pPi=REAL(pi);
+ int *pA1=INTEGER(A1), *pA2=INTEGER(A2);
+ double *pRho=REAL(rho), *pPi=REAL(pi);
/////////////////////////////
// Call to selectiontotale //
/////////////////////////////
- selectiontotale(
+ selectiontotale_core(
phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,
pA1,pA2,pRho,pPi,
n,p,m,k,L);
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);
+ SET_VECTOR_ELT(listParams, 0, A1);
+ SET_VECTOR_ELT(listParams, 1, A2);
+ SET_VECTOR_ELT(listParams, 2, rho);
+ SET_VECTOR_ELT(listParams, 3, pi);
UNPROTECT(7);
return listParams;
#include <gsl/gsl_linalg.h>
// TODO: don't recompute indexes every time......
-void EMGLLF(
+void EMGLLF_core(
// IN parameters
const double* phiInit, // parametre initial de moyenne renormalisé
const double* rhoInit, // parametre initial de variance renormalisé
#ifndef valse_EMGLLF_H
#define valse_EMGLLF_H
-void EMGLLF(
+void EMGLLF_core(
// IN parameters
const double* phiInit,
const double* rhoInit,
}
// TODO: comment EMGrank purpose
-void EMGrank(
+void EMGrank_core(
// IN parameters
const double* Pi, // parametre de proportion
const double* Rho, // parametre initial de variance renormalisé
#ifndef valse_EMGrank_H
#define valse_EMGrank_H
-void EMGrank(
+void EMGrank_core(
// IN parameters
const double* Pi,
const double* Rho,
#include <omp.h>
// TODO: comment on constructionModelesLassoMLE purpose
-void constructionModelesLassoMLE(
+void constructionModelesLassoMLE_core(
// IN parameters
const double* phiInit, // parametre initial de moyenne renormalisé
const double* rhoInit, // parametre initial de variance renormalisé
double* piLambda = (double*)malloc(k*sizeof(double));
double* LLF = (double*)malloc((maxi+1)*sizeof(double));
double* S = (double*)malloc(lengthA*m*k*sizeof(double));
- EMGLLF(phia,rhoInit,piInit,gamInit,mini,maxi,gamma,0.0,Xa,Y,tau,
+ EMGLLF_core(phia,rhoInit,piInit,gamInit,mini,maxi,gamma,0.0,Xa,Y,tau,
phiLambda,rhoLambda,piLambda,LLF,S,
n,lengthA,m,k);
free(Xa);
for (int mm=0; mm<lengthB; mm++)
{
for (int r=0; r<k; r++)
- phi[ai( A2[ai4(j,0,lambdaIndex,p,m+1,L)]-1, b[mm], r, lambdaIndex, p, m, k, L)] = 0.0;
+ phi[ai4( A2[ai(j,0,lambdaIndex,p,m+1,L)]-1, b[mm], r, lambdaIndex, p, m, k, L)] = 0.;
}
}
#ifndef valse_constructionModelesLassoMLE_H
#define valse_constructionModelesLassoMLE_H
-void constructionModelesLassoMLE(
+void constructionModelesLassoMLE_core(
// IN parameters
const double* phiInit,
const double* rhoInit,
#include "utils.h"
// TODO: comment on constructionModelesLassoRank purpose
-void constructionModelesLassoRank(
+void constructionModelesLassoRank_core(
// IN parameters
const double* Pi,// parametre initial des proportions
const double* Rho, // parametre initial de variance renormalisé
for (int v=0; v<m; v++)
{
for (int r=0; r<k; r++)
- RhoLambda[ai(uu,v,r,m,m,k)] = Rho[ai4(u,v,r,lambdaIndex,m,m,k,L)];
+ RhoLambda[ai(u,v,r,m,m,k)] = Rho[ai4(u,v,r,lambdaIndex,m,m,k,L)];
}
}
- EMGrank(PiLambda,RhoLambda,mini,maxi,Xactive,Y,tau,rank,
+ EMGrank_core(PiLambda,RhoLambda,mini,maxi,Xactive,Y,tau,rank,
phiLambda,&LLF,
n,longueurActive,m,k);
free(rank);
#define valse_constructionModelesLassoRank_H
// Main job on raw inputs (after transformation from mxArray)
-void constructionModelesLassoRank(
+void constructionModelesLassoRank_core(
// IN parameters
const double* Pi,
const double* Rho,
#include "utils.h"
// Main job on raw inputs (after transformation from mxArray)
-void selectiontotale(
+void selectiontotale_core(
// IN parameters
const double* phiInit, // parametre initial de moyenne renormalisé
const double* rhoInit, // parametre initial de variance renormalisé
#define valse_selectiontotale_H
// Main job on raw inputs (after transformation from mxArray)
-void selectiontotale(
+void selectiontotale_core(
// IN parameters
const double* phiInit,
const double* rhoInit,
projects / valse.git / commitdiff
